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The increase in fuel prices has initiated considerable 
interest by ship operators in new ship autopilots which 
minimize the propulsion losses due to steering. 

This research presents the results of work on a steering 
design for the Mariner class ship based on a computer simu- 
lation. A model of the Mariner class ship was coupled to a 
function minimization subroutine to minimize the added 
resistance caused by rudder activity and hull drag of iner- 
tial origins caused by periodic yawing of ship in seaway. 

The Mariner class ship computer model was tested in calm 
water and in a seaway. The optimal controller parameters 
are shown in look up tables as functions of ship speed, sea 
state, encounter angle and encounter frequency. This tech- 
nique can be used as' an adaptive controller. 
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I . INTRODUCTION 



Claims by many researchers indicate that a carefully 
designed controller could reduce fuel consumption by mini- 
mizing the propulsion losses which are caused by added drag 
due to steering of the ship. 

The goal of this thesis was aimed at developing and 
demonstrating the utility of an improved steering control 
for the Mariner class ship. The immediate goal of this 
research was to develop design methodology for an adaptive 
autopilot that would provide effective steering control with 
associated cost savings for a full range of seaway and 
stability conditions. 

To simulate the ship in the computer program, ship's 
nonlinear equations of motion were needed. Chapter 2 
addresses the Mariner class ship nonlinear equations of 
motion. 

Chapter 3 addresses the work 
tion model and open loop ship's 
in a seaway. 

The basic Nomoto models give an adequate description of 
ship steering dynamics for design. The Nomoto second and 
third order models were developed from the ship's linear 
equations of motion, Chapter 4 adresses the Mariner class 
ship Nomoto model derivations by using mathematical methods 
and by using a function minimization subroutine. 

Chapter 5 shows an adequate cost function which repre- 
sents the added drag due to steering and includes deriva- 
tions for evaluating the weighting factor. Also Chapter 5 



on testing the ship simula- 
behaviors in calm water and 
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presents the assumptions and approaches needed to find the 
cost function that is used by many researchers. 

Ship dynamics change with operating conditions such as 
ship speed, encounter angle and encounter frequency. 
Chapter 6 presents optimal controller parameters as a func- 
tion of different operating conditions. 

Chapter 7 addresses an approach to an adaptive controller 
utilizing information which is easy to measure on ship board 
such as ship speed, heading error and rudder angle. This 
adaptive controller must be used to provide minimum added 
drag due to steering. 

Conclusions were drawn from simulation results and are 
presented in Chapter 8. This chapter also recommends some 
future studies, which can be done as extensions of this 
work. 
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II. NONLINEAR EQUATIONS OF MOTION 



Nonlinear equations of motion are suitable for predicting 
tight maneuvers and also suitable for computer programming. 
The nonlinear equations of motion used in this work have 
been developed by Abkowitz [Ref. 1, 2], and Strom_Tejsen 

[Ref. 3], based on a Taylor series expansion of forces and 
moments. Terms higher than third order are not included in 
the equations because experience has shown that accuracy is 
not significantly improved by their inclusion. 

A result of symmetry about the xz-plane, X is an even 

» • 

function of V, R, D, V and R so on, the crosscoupled terms 
in the equations involving odd powers of V, R and D are 
zero, however, crosscouple terms which involve even powers 

of V,R and D are nonzero. In contrast to X, the expressions 

« • 

for Y and N are odd functions of V,R,D,V and R; that is, 
only the coefficients of the terms in the expansion with odd 
powers are nonzero; those with even powers are zero. For 
some reasons , X is neither an odd nor an even function of U 
but rather its expansion includes all powers of DU. 

• • • 

Equations X,Y and N are functions of U,V,R,U,V,R and D. 
Taylor series expansions of X,Y and N including terms up to 
the third order are as follows: 



(m - X^)*U = f/U,V,R,D) 



(eqn 2.1) 



(m - Y^)*V + (m*X^- Y^)^^R = f/U,V,R,D) 



(eqn 2.2) 



(m*X^- N. )^^^V + (I^- N-)*R = f^(U,V,R,D) 



( eqn 2.3) 
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Where : 



d^X/dU^ X^^= d'X/dD^ X^^= d^X/dU^’^dV 

X^^^ d’X/dU^ X^^= d’X/dU^*dR etc. 

f^(U,V,R,D) = 0.5*X^^'^DU^ + (1/6)*X^^^*DU' (eqn 2.4) 

+ 0.5*Xyv*V^ + (0.5*X^^+ m*X^)*R" 

+ 0.5*Xyvu;-V^*DU + 0.5“X^^.^-R"“DU 
+ (X^^+ m)*V*R + X^o*V*D + X^i^^'^V^R^DU 
+ Xvova*V“D"DU + X^^^'-R*D*DU + X^-DU + X“ 

+ + 0.5^^X,^^'D^*DU ^ X^^^^R^^D 



f^(U,V,R,D) 



= Y“^ -VDU + Y^^*DU’^+ Yv'*V + Yq*D (eqn 2.5) 

+ Y“ + 0.5*Y^p_^*V*R^ + 0.5 “Y^ot>'‘V"D" 

+ Yvvj:’'V"DU + 0.5*Y^v,,vj:’'V-'DU^ + (Y^- m*Ul)*R 
+ (1/6)*Y,^^ + Yr.^^-^R^^DU + 0.5*Y^^^^R*DU^ 

+ (1/6)“Y^^^*R’ + + 0 . 5*Y^yv*R*V 

+ 0.5*Yjj(^p^*D“R^ + Yo^*D*DU + 0 . 5*Y^^^*D'--DU^ 

+ Y^(^o*V'^R*D + (1/6)*Y^^^*V’ + 0.5 *Y^jjd*R"D^ 



2 



f(U,V,R,D) = “DU + Nt,vx*DU^ + N^*V + Np *D (eqn 2.6) 



+ 0.5*Np^^*D^^DU^ 



(Np-m^'^X(iUl)*R 



+ N^^*V*DU + 0.5*N^^vj:-V“DU" + N^p_q*V*R"D + N" 
+ 0. 5*Nq^^*D*V" + 0.5“Nt)|^R^*D“R^ + N^^-'D*DU 
+ N^^-^R^^DU + O-S^N^^^^R^-^DU^ + ( 1/ 6 )*N^jj/D ^ 

+ (1/6)*N^^^*R' + + 0 . 

+ (1/6)^^N^^/V’ + 0.5*N^,^p;’^V^'^R^ + 0 . 



All of the derivative coefficients of the equations are 
evaluated on the basis of experimental data obtained from 
captive model tests, and given in Table I, II and III. 

The Y force and N moment induced by the rotation of a 
single propeller or by unirotating multiple propellers at 
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0 o 

V=D=0.0 are identified as Y and N, these are likely to be 
speed dependent. To see the rudder effects in calm water or 
sea state propeller effects were ignored. 

Finally, from the X , Y and N equations the ship's surge, 
sway and yaw equations can be written as follows. 



f^(U,V,R,D) 



U = 



(eqn 2.7) 



(m - X^) 



(I^- Nj^)^-^f^(U,V,R,D) 



(m^-^X^- Y^)^^f^(U,V,R,D) 



V = 



( eqn 2.8) 



(m - Y.)*(I^- N^) - (m^'^X^- Y- ) 



(m - Y;^)'-^f^(U,V,R,D) 



(m^-^X^^- N- )*yU,V,R,D) 



R = 



(eqn 2.9) 



(m - Y. )*(I^- N^) - (m*X^- Y^) 
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These surge, sway and yaw equations can be rewritten in the 
f orm : 

dU/dt = g^(t, U(t), V(t), R(t), D(t)) 

dV/dt = g^(t, U(t), V(t), R(t), D(t)) (eqn 2.10) 

dR/dt = g^(t; U(t), V(t), R(t), D(t)) 

Where U(t), V(t), R(t) and D(t) are the instantaneous values 
of U, V, R and D at any time t. 

Equation 2.10 is a set of three first-order differential 
equations for which approximate numerical solutions are 
readily obtained on a digital computer. The key to the 
numerical solution is that values of U, V and R at time 
t+DELT are obtained from knowledge of the values of U, V, R 
and D at time t using a simple first-order expansion; that 
is , 

U(t+DELT) = U(t) + DELT“U(t) 

V(t+DELT) = V(t) + DELT"V(t) (eqn 2.11) 

R(t+DELT) = R(t) + DELT*R(t) 

This method is found to give adequate accuracy for the 

present type of differential equations because of the fact 

• • • 

that the accelerations U, V and R vary but slowly with time, 
due to the large mass and inertia of a ship compared to the 
relatively small forces and moments produced by its control 
surface. Any desired accuracy of the solutions can be 
obtained with a computer by simply using smaller time inter- 
vals DELT. This procedure was used for all computer programs 
which were developed for this thesis. 
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TABLE I 

Assessment of the Coefficients^ in the X_Equation 



Variable 


Coefficient 


Value of Coeff: 


U 


m - 


0.177 


DU 




-0.0253 


DU^ 


0.5"X^^ 


0.00948 


DU^ 


(1/6)"X^^^ 


-0.00217 




0.5-X^,, 


-0.189 


R' 




0.00379 






-0.02 


V^*DU 


0 . 5 





R^"DU 


0 . 5 "X 


— 


D"*DU 


0.5-^X^^^ 





V"R 


XvR.^ m 


0.168 


V"D 


Xs,o 


0.0196 


R"D 


XRt> 


0 


V*R*DU 


X Vgui. 





V-D"DU 


X V Ov*. 





R*D--DU 


X 







X® 


0 



^All derivatives are nondimensionalized on the basis of 
RHO,L,T and S. 

No entry in these columns means the coefficient was ignored. 
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TABLE II 

Assessment of the Coefficients in the Y_Equation 



Variable 


Coefficient 


Value of Coefficient 


V 


(m - Y^) 


0.327 


• 

R 


Y^) 


-0.0018 


V 




-0.244 


V’ 


(1/6)*Y^^^ 


-1.702 


V*R2 




0 


V"D^ 




-0.0008 


V-DU 




_ 


V"DU^ 


0 5*Y 





R 


(Y^- m) 


-0.105 


R^ 




0 


R " V ^ 


0-5^^Y^vv 


3.23 


R-D^ 




0 


R*DU 


Y^vx 


_ 


R”DU^ 







D 


Yd 


0.0586 


D’ 


(1/6)’^Y^^^ 


-0.00975 




0.5'-'Y^^^ 


0.25 


CM 

Q 


0 . 5"Y^^^ 


0 


D-'DU 


Ytvx 


0 


D " DU ^ 






Q 

<* 

> 


Yn/RO 


0 





Y“ 


-0.0008 


DU 


Y“ 


0 


DU^ 


yO 

^ \xvjv 
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TABLE III 

Assessment of the Coefficients in the N_Equation 



Variable 



Coefficient 



Value of Coefficient 



V 


(m*X^- N-) 


-0.00478 


% 

R 


(I^- N^) 


0.0175 


V 


Nv 


-0.0555 


V^ 




0.345 


V'-R" 


0 . 


0 


V-'D^ 




0.00264 


V*DU 


Nvvx 


_ 


V*DU ^ 


0.5*Nvu.si, 


_ 


R 


(N^- m^^X^-Ul) 


-0.0349 


R' 


(1/6)^^N^^^ 


0 


R"V" 




-1.158 


R-D" 




0 


R"DU 


Ni^vx 





R “ DU " 







D 




-0.0293 


D' 




0.00482 


D*V" 


0 . 5 


0.1032 


D*R" 




0 


D"DU 




0 


D*DU^ 


0-5-Nw^^ 


— 


V*R"D 


^'JRT) 


0 




N“ 


0.00059 


DU 


N" 


0 


DU" 


N° 

\x\x 
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III. OPEN LOOP SHIP'S BEHAVIOR IN CALM WATER AND SOME SEA 

STATES 



A. CALM WATER CASE 

Using a ship's nonlinear equations of motion, and Mariner 
class ship coef f icients , a simulation program THESIS FORTRAN 
was developed and run to observe U, V and R.The computer 
program THESIS FORTRAN is shown in Appendix A. 

First run D=0, YAWC=0, Ul=15 knots were applied and it 
was seen that the ship stays on its initial course and 
speed. U=15 knots, V=R=0. 

The program was rerun a few times, changing the rudder 
angle to 2.8 degrees to both sides (port and starboard) and 
it was observed that increasing the rudder angle changes the 
ship's course and also U decreased, the absolute values of V 
and R increased. After a few hundred seconds U, V and R 
reached steady values independently. These steady values 
depend on rudder angle, large rudder angles decrease U and 
increase V and R. For a rudder angle of one degree and 
constant speed (Ul) of 15 knots, the first 200 seconds of 
time response of U, V and R are shown in figure 3.1 as an 
example . 

B. SEA STATE CASE 

To observe the behavior of the ship in a sea state, 
disturbance forces and moments are needed that depend on sea 
state, ship speed, encounter angle and encounter frequency. 
Also in sea steate hydrodynamic parameters are changed, i.e, 
the added mass and added inertia are functions of encounter 
frequency and sea state. 
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The disturbance forces, moments, and added mass, added 
inertia terms were found by running the sea state program 
that is presented in [Ref. 4]. Data for a Mariner class 
ship and chosen sea conditions that were used in the sea 
state program are shown in APPENDIX B. 

For observation purposes ' FX ' (disturbance force in 
surge), 'FY* (disturbance force in sway) and 'MZ' (distur- 
bance moment in yaw) were added into the surge, sway and yaw 
equations that were used in the simulation program. 

For regular seas the program has been run a few 
times. Every time different FX, FY, MZ and coefficients were 
used to represent different ship characteristics such as 
ship speed, loading ^ etc. and enviromentel conditions such 
as sea state, encounter angle, encounter frequency. Results 
show that U, V and R are sine waves with amplitude and phase 
depending on ship and environmental conditions. 

Time response of U, V and R in 200 seconds are presented 
in Figure ( 3.2), for ship speed 15 knots, encounter angle 
030.0 degree, encounter frequency 0.60 radian/ second and sea 
state 6. 



^During this research, displacements (Molded) up to 32 
ft were chosen as a loading condition for the Mariner. 




O o 



TIME (SEC.) 




Figure 3.1 Time Response of U, V AND R in Calm Water when 

D = 1 Degree . 
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TIME (SEC.) 




TIME (SEC.) 




Figure 3.2 Time Response of U,V and R in Regular Seas. 
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IV. NOMOTO MODEL OF THE MARINER 



To find a model which can be used in computer simulation 
the Mariner's linear equations of motion and its hydrody- 
namic coefficients were used and third and second order 
Nomoto transfer functions derived. Values used were for ship 
speed of 15 knots. 

Mariner's linear equations of motion are 

(m - = X^*DU (eqn 4.1) 

(m - Y-y-'V - Y^*V = (Y— m*Xc)*R + (Y,^- m*Ul)*R 
(I - Np^4 - (Nj^- m*X^*Ul)*R = (N<,- m*X^^)*V + N^*V 

A. MATHEMATICAL APPROACH 

Proceeding to the third order Nomoto equation: 

YAW(s) K*(Z-s + 1) 



D(s) s*(Pl*s + l)*(P2*s + 1) 

Deriving the third order Nomoto transfer function from 
the sway and yaw equations, we show them in matrix notation 
as follows. 









— 




— 




— — 


V 




- .0372 


-8.42 




V 




.210 




r 






* 








R 




- .0003 


- . 10 




R 




- .003 






— 


_ 


, 






— 
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X = A*X + B"U then, X(s) 
conditions are zero.) 



( s*I- A)"^ (assuming initial 



V 




.2101*s + 


.0534 


R 




-.0038*s - 


.0002 



D 



(7.67*s+l)*(116.93*s+l) 



R(s) = s*YAW(s) 



YAW(s) ■ -0.189*(18.34*s+l) 



so, then 



D ( s ) s* ( 7 . 6 7 - s + 1 ) * ( 116 . 9 3*s + 1 ) 



result is K=0.189, z=18.34, Pl=7.67 and P2=116.93 
Proceeding to the second order Nomoto equation: 
YAW(s) K 



D(s) s''(Pl*s + 1) 

Deriving the second order Nomoto transfer function from 
the yaw equation only, the result is K=0.03 and Pl=10. 

B. COMPUTER APPROACH 

We used a function minimization subroutine to obtain 
parameters of the transfer functions. Figure 4.1 shows the 
scheme used to obtain the third and second order Nomoto 
transfer functions. The results are given in Table IV. 
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+ 




Figure 4.1 Determination of Third and Second Order 

Nomoto Models . 
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TABLE IV 

The Nomoto Model Parameters for Mariner 



Third Order 


Second Order 


K = 0.189 


K =0.0298 


Z = 18.347 
PI = 7.6739 
P2 = 116.929 


PI = 9.989 



As seen the answers obtained by function minimization 
agree closely with the analytic solutions. 
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V. DERIVATION OF A COST FUNCTION FOR THE MARINER SHIP 



In recent years many researchers have studied the problem 
of optimizing an automatic ship steering controller for 
minimum fuel consumption. It is well known that additional 
drag is introduced by steering and that both the rudder 
motion and the yawing motion contribute to this added drag. 

When deriving a cost function we are required to find one 
cost function that must be convenient for ship board use. The 
cost function that is commonly used in recent years is 



Derivation of this cost function from the surge equation has 
been well explained in [Ref. 5] for the SL/7 class ship by 
R.E.REID. 

To derive the cost function for the Mariner, REID's 
approach is taken as a reference and his assumptions are 
used. To show how the cost function for the Mariner was 
derived, REID's work is presented here step by step with 
derivations of the Mariner's cost function. 

The Surge Equations: 




( eqn 5,1) 



0 



Where YAWE = Yaw error 

LAMBDA = Weighting factor 



SL/7 : 




+ (Xyg_ + m)“V“R + Xp 
Where: X® = -0.0003 for SL/7 



Xp = Propeller thrust 
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MARINER; 



(m - X^)^-U = 0.5*X^^*V^ + (XvR. ^ m)’^V“R (eqn 5.3) 

+ 0 . 5 *X '• D 2 + (0.5 -X + m " X ^ ) " R * 

+ Xvo"V"D + X^“DU + 0.5''X^^'^DU^ 

+ (1/6)^^X^^^^DU’ 



It is seen that there are some terms in the Mariner's 
surge equation which they are not included in the SL/7's 
surge equation. Assuming steady state situations since, U=0. 
The instantaneous surge relevant to steering is 



DX = 0.5''X^^*V" + 0.5-’'X^^''D" + (X^^+ m)*V*R (eqn 5.4) 



SL/7 : 

MARINER: 

DX = 0.5“Xvm"V" + 0.5-'Xx)X)*D" + + m)“V"R (eqn 5.5) 

+ (0.5-X^^ + m*X )''R" + Xs,t,'"V"D + X^’'DU 



From the instantaneous surge equation relevant to 
steering of the SL/7 Reid came up with the following cost 
function: 



J = (1/2T)* / (LAMBDA "*V"R + ETA*V" + D^'dt (eqn 5.6) 



0 

Where 



LAMBDA"= (m + X^^) / 0 . 5^*^X^^ 
ETA = (0.5*X,^,,)/0.5*X^.5^ 
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and he found that the value of the (ETA'-V^) term is very 
small so he neglected the (ETA"V*) term, then the cost func- 
tion for the SL/7 is 



•t.. 



J = 0. 5* /(LAMBDA" *V*R + D*)*dt 



(eqn 5.7) 



Using the same approach for the Mariner, the following 
cost function was derived. 

J = (1/2T)*/(A1-DU + A2"DU" + A3 -'DU' + A4*V" (eqn 5.8) 



A5"R" + + 


A7 


“V*R - A8*V*D)*dt 






Where Al 


= 


X^/0.5^^X^^ 


A2 = 


0.5*X^,,/0.5*X^^ 


A3 


= 


(1/6)*X^^^/0.5*X^^ 


A4 = 


0.5*X^^/0.5*X^^j 


A5 


= 


(0.5-X^^ + m"Xg)/0. 






A7 


= 


(X,^^ * m)/0.5*X^^ 


A8 = 


X^o/0-5“X^;, 



A5 and A8 are always negative numbers, since then A5 and A8 
have a minus sign in the equation. For the calm water case 
and when Ul=15 knots, D=2.6 degrees after 2000 seconds, 
values of every term in the Mariner's cost function are 
given below to give an idea about assumptions. 



Al-DU = 


-0.001939 


A2‘-DU^ = 


-0.00000111 


A3-DU*= 


-0.00000000039 


A4*V^ = 


0.000003253 


A5*r2 = 


-0.0002884 


A7*V'>R = 


0.0001923 




0.002059 


A8*V“D = 


-0.00002609 



As seen from the above (A2*DU^), (A3“DU^), (A4*V^ ) and 

( A8"'V”D) terms are very small compared to others, so they 
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may be neglected. Also to measure the DU on shipboard is 
very hard although it may someday be measured by new satel- 
lite facilities, so we must not include it in the cost func- 
tion. After these assumptions the cost function for the 
Mariner is 



A. 



J = 0.5" (-A5"R^ + D" + A7"V"R)"dt 



( eqn 5.9) 



0 

The only difference between Equation (5.9) and Equation 
(5.7) is the (A5"R^ ) term that is not included in the SL/7's 
cost function. To see the effect of the (A5"R^ ) term on the 
cost function the Mariner's cost function was evaluated with 
and without the (A5"R^ ) term for the calm water case and 
Ul=15 knots, D=2.6 degree after 2000 seconds, results are: 



J 



with (A5"R^ ) 


0.0025399826296 


without (A5"R^ ) 


0.0022515066462 



There is no big difference between these two J values, and 
to make the derivations similar to Reid's derivations the 
(A5"R^ ) term won't be included in the Mariner's cost func- 
tion but as it is known that for the Mariner the (A5"R^ ) 
term is as big as the (A7*V"R) term, it would be better to 
consider it in the cost function. After all of the above 
steps the cost function of the Mariner may be written as in 
Equation (5.7). 
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V and R are hard to measure on shipboard, but (V'-R) can 
be defined as 

V*R = OP -YAW ^ 

Where R = YAW = YAWE-w 

OP = Distance from the ship pivot point to 
the origin. = 0.3*L 

w = Natural frequency of the ship's steering 
system closed loop, (w = 0.05 rad/sec was initially used.) 

Finally the cost function for the Mariner is 



J = 0 . 5 •' / ( LAMBDA -YAWE ^ + D ^ ) *d t 



( eqn 5.10) 



Where LAMBDA = (m + )*0P'''w VO . 5*X 






Since, X depends on ship speed, for different ship speeds 
values of LAMBDA were calculated and are presented in Table 
(V). These LAMBDA values were calculated by assuming the 
natural frequency of the ship's steering system closed loop 
is equal to 0.05 radian/ second . How important the accuracy 
of the LAMBDA value is with respect to finding the optimal 
control parameters will be observed in Chapter 6. 
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TABLE V 



Values of Weighting Factor for Different Speeds 
of The Mariner Class Ship 



ship speed | 


LAMBDA 


(Knots) 1 




10 1 


6.57 


15 1 


2.91 


20 1 


1.64 
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VI. CONTROLLER DESIGN FOR THE MARINER. USING FORTRAN PROGRAM 



We coupled a function minimization subroutine to the 
cascade compensater that is coupled with Mariner’s equations 
of motion and used the subroutine to adjust the controller 
parameters to minimize the cost function which is derived in 
Chapter V, and evaluate the minimum cost. The Fortran program 
to calculate the optimal parameters is given in Appendix C. 
Compensater 'A' and 'B' are used as controllers, their 
structures are shown in Figure ( 6.1 ). 



Kx(SxZl f 1) 




Kx(SxZl + 1) 


(SxPI f 1 ) 




(SxPI f 1)x(SxP2 f 1) 



COMPf N'JAIOR - A COMPENSATOR 8 



Figure 6.1 Various Structure Controllers. 



A. CALM WATER CASE 

For calm water, for a given 1 degree yaw command, using 
the computer method we optimized controllers ’A' and 'B' and 
the results are shown in Tables VI, VII and VIII 
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TABLE VI 

Optimal Controller Parameters I 



Simulation Results - Steady State 600 seconds 

For ship speed 10 Knots the optimal 
parameters of various controllers 
and the cost 



CONTR 


K 




PI 


P2 


COST 


A 


0.77 


60.07 


20.46 




0.0494896 


B 


0.74 


60.07 


20.12 


0.902 


0.0514841 



TABLE VII 

Optimal Controller Parameters II 



Simulation Results - Steady State 600 seconds 

For ship speed 15 Knots the optimal 
parameters of various controllers 



CONTR 


and the cost 
K Z1 


11 


?2 


COST 


A 


0.53 


51.46 


18.08 




0.0186974 


B 


0.50 


51.58 


17.81 


0.890 


0.01957901 



TABLE VIII 

Optimal Controller Parameters III 



Simulation Results - Steady State 600 seconds 



For ship speed 20 Knots the optimal 
parameters of various controllers 
and the cost 



CONTR 


K 


?1 


PI 


P2 


COST 


A 


0.40 


44.87 


16.06 




0.0094217 


B 


0.39 


41.11 


15.84 


0.880 


0.00991801 
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These results will be references for the controller 



design for sea state operation. We observe that increasing 
the speed gives us smaller controller parameters, also this 
behavior can be seen from the SL/7's results that are 
presented in [Ref. 6]. 

B. REGULAR SEAS CASE 

The ship in regular seas is affected by sea wave distur- 
bance forces and moments. These are functions of sea state 
and encounter frequency. Also the added mass and the added 
inertia terms are functions of sea state and encounter 
frequency, and the encounter frequency depends on the 
encounter angle and ship speed. All of these variables must 
be considered when calculating the optimal parameters of the 
controller . 

Using the computer method controllers 'A' and 'B' were 
optimized for a few different cases and the results are 
shown in Tables IX, X, XI, XII, XIII, XIV, XV, XVI and XVII. 

TABLE IX 

Optimal Controller Parameters IV 

Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 

For ship speed 15 Knots, encounter angle 030.0 degree, 
encounter frequency 0.50 rad. /sec. and sea state 6. 



CONTR 


K 




£1 


P2 


COST 


A 


0.358 


66.6 


24.61 




0.001252939 


B 


0.35 


44.68 


06.58 


9.720 


0.0010086581 
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TABLE X 

Optimal Controller Parameters V 



Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 



For ship speed 15 Knots, encounter angle 060.0 degree, 
encounter frequency 2.50 rad. /sec. and sea state 9. 



CONTR 


K 




PI 


P2 COST 


A 


0.60 


41.58 


11.33 


0.000020747 


B 


0.70 


30.49 


03.37 


2.940 0.000016038 



TABLE XI 

Optimal controller Parameters VI 



Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 

For ship speed 15 Knots, encounter angle 090.0 degree, 
encounter frequency 0.50 rad. /sec. and sea state 7. 



CONTR K ^ PI P2 COST 

A 0.54 30.65 09.35 0.054223988 



g XT DID NOT CONVERGE 



TABLE XII 

Optimal Controller Parameters VII 



Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 

For ship speed 15 Knots, encounter angle 120.0 degree, 
encounter frequency 0.75 rad. /sec. and sea state 7. 



CONTR K ^ PI ^ COST 

A 0.52 41.09 08.95 0.018345 



B *.v*****.v** XT DID NOT CONVERGE 
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TABLE XIII 

Optimal Controller Parameters VIII 



Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 

For ship speed 15 Knots, encounter angle 150.0 degree, 
encounter frequency 1.50 rad. /sec. ana sea state 8. 



CONTR K 



Z1 PI P2 COST 



A 

B 



0.645 




42.40 07.70 

IT DID NOT CONVERGE 



0.0008188 



TABLE XIV 

Optimal Controller Parameters IX 



Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 

For ship speed 15 Knots, encounter angle 090.0 degree, 
encounter frequency 0.50 rad. /sec. and sea state 8. 



CONTR 



K 



Z1 



PI 



P2 



COST 



A 

B 



0.53 



27.45 08.04 

IT DID NOT CONVERGE 



0.070697939 



TABLE XV 

Optimal Controller Parameters X 

Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 

For ship speed 15 Knots, encounter angle 090.0 degree, 
encounter frequency 0.50 rad. /sec. and sea state 6. 



CONTR K ^ Pf £2 COST 

A 0.56 39.71 11.87 0.019029424 



B XT did NOT CONVERGE 
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TABLE XVI 

Optimal Controller Parameters XI 

Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 



For ship speed 10 Knots, encounter angle 060.0 degree, 
encounter frequency 2.50 rad. /sec. and sea state 9. 



CONTR 


K 




Zi 


11 


COST 


A 


0.80 


45.00 


10.66 




0.000045194 


B 


0.89 


36.27 


03.22 


03.22 


0.000036077 



TABLE XVII 

Optimal Controller Parameters XII 



Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 



For ship speed 20 Knots, encounter angle 060.0 degree, 
encounter frequency 2.50 rad. /sec. and sea state 9. 



CONTR 


K 


zi 


PI 


il 


COST 


A 


0.47 


34.67 


10.35 




0.000014376 


B 


0.61 


23 . 16 


03.23 


03.03 


0.000011085 



As seen from the above results, for all cases the 
compensators have characteristics of a lead network. 

The effects of sea state on the controller parameters 
can be seen by comparing the Tables XV, XI and XIV As seen 
from the tables an increase in sea state causes an increase 
in the cost value and a decrease in the controller parame- 
ters as expected, because heavy sea state brings high 
disturbance forces and moments and they cause heavy yawing 
motions. To answer this, the controller time constants 
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decrease and new parameters introduce more rudder motion, so 
an increase in rudder motions and heading error causes 
increase in the cost. For ship speed 15 knots, encounter 
angle 90.0 degree, encounter frequency 0.5 radian/ second , 
sea state 6 and sea state 8, rudder angles and heading 
errors were plotted in 200 seconds and they are shown in 
Figure ( 6.2 ) and ( 6.3). 
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o 




TIME (SEC.) 



Figure 6.3 Heading Errors in Degrees for Sea State 6 and 8. 

It can be seen from Figure ( 6.2 ) and ( 6.3 ), for sea 
state 8 rudder and heading error values are bigger than the 
sea state 6 rudder and heading error values, but it was 
noticed that the periods are almost the same for both sea 
states . 
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To have a better understanding about the effects of sea 
state on optimal controller parameters the cost value versus 
parameter values were plotted for ship speed 15 knots, 
encounter angle 90.0 degree, encounter frequency 0.50 
radian/ second and sea state 6 between 8 and they are shown 
in Figure ( 6.4 ), ( 6.5), ( 6.6) and ( 6.7). 




K 



Figure 6.4 The Cost vs. K. 
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0.01 0.03 0.05 0.07 




Z1 



Figure 6.5 The Cost vs. Zl. 
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PI 



Figure 6.6 The Cost vs. PI. 
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Figure 6.7 The Cost vs. Z1 and PI. 



As seen from the figures the optimal parameter values 
have a linear behavior when changing the sea state. This 
property would be useful to create look up tables for param- 
eter values. The use of look up tables will be discussed in 
Chapter 7. 

To see the motion of the ship, the ship's equations are 
solved in ship coordinates (x , y) and then transformed to 
space coordinates (x , y), Figure (6.8 ) shows the orienta- 
tion of space axes and moving axes. 
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Figure 6.8 Orientation of Space Axes and Moving Axes. 



The transformation from ship to space coordinates is 
defined as 



Xog = U*cos(YAW) - V" sin (YAW) (eqn 6.1) 

Yog = U“ sin (YAW) - V" cos (YAW) 



Xog, Yog = Coordinates of the center of mass of 

the ship relative to coordinate system 
fixed with respect to the surface of the 
earth 



A 6 



Where 



Equation ( 6.1 ) was added into the fortran program and 
the ship's motion were observed for 600 seconds in calm 
water and in regular seas for 30.0 degree encounter angle, 
sea state 6, encounter frequency 0.60 radian/ second and ship 
speed 15 knots. Both compensator 'A' and 'B' results are 
shown in Figure ( 6.9 ), and after 600 seconds ship's coor- 
dinates are given in Table XVIII. 



TABLE XVIII 
Location of the Ship 

Simulation Results - Steady State 600 seconds 

CASE Xog (ft) Yog (^) 

Calm water 

No rudder 15200.994873 0.00 



Regular seas 

Comp. 'A' 15199.0237114 3.2596 

Regular seas 

Comp.'B' 15199.0763219 2.87116 



For calm water and regular seas compensator 'B' gave 
better results than compensator 'A' did, but for some cases 
compensator 'B' did not converge and the difference between 
results is not significant. The comparisons were made as to 
which type of compensator brings the ship nearest to final 
location at the end of the 600 seconds. 
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YOG (K'l') 

2.5 




Figure 6.9 Simulation Results- Steady State 600 sec. 
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To see the effects of the weighting factor ( LAMBDA ) , 
on the compensator parameters, different LAMBDA values, were 
used with compensator ’A’ and the compensator parameters 
were computed for ship speed 15 knots, calm water and 
regular seas with encounter angle 90.0 degree, encounter 
frequency 0.5 radian/ second and sea state 8. Results are 
shown in Table ( XIX) and Table ( XX ) . 



TABLE XIX 

Test for LAMBDA Value I 

Simulation Results - Steady State 600 seconds 

for different weighting factor values , 
calm water, ship speed 15 knots, 
optimal parameters of A type controller. 



LAMBDA 


K 




PI 


10.00 


1.07 


29 . 59 


11.30 


9.00 


1.00 


31.5 


11.90 


8.00 


0.94 


32.91 


12.45 


7.00 


0.88 


34.84 


13.14 


6.00 


0.80 


37 .61 


14.02 


5.41 


0.76 


39.38 


14.62 


4.91 


0.72 


41.26 


15.19 


4.41 


0.68 


43.14 


15.78 


3.91 


0.63 


45.26 


16.36 


3.41 


0.58 


48.06 


17 . 10 


2.91 


0.53 


51.07 


17 . 93 


2.41 


0.48 


55.21 


18 . 99 


1.91 


0.42 


59 . 79 


20.15 


1.41 


0.35 


64.79 


21.31 


0.91 


0.27 


76.89 


24.17 


0.41 


0.18 


89 . 61 


28.22 
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TABLE XX 

Test for LAMBDA Value II 



Simulation Results - Steady State 600 seconds 

for different weighting factor values, 
regular seas, ship speed 15 knots 
encounter angle 90.0 degree, sea state 8, 
and encounter frequency 0.50 radian/ second 
optimal parameters of A type controller. 



LAMBDA 


K 


21 


£1 


5.41 


0.71 


24.91 


08.41 


4.91 


0.68 


25.36 


08.37 


4.41 


0.64 


25.88 


08 . 32 


3.91 


0.60 


26.34 


08.22 


3.41 


0.56 


26.85 


08.11 


2.91 


0.53 


27.45 


08.04 


2.41 


0.48 


28.28 


08.04 


1.91 


0.44 


29 . 16 


08 . 16 


1.41 


0.40 


31.21 


08.35 


0.91 


0.34 


33.95 


08 . 71 


0.41 


0.26 


39.00 


09.27 



Using above compensator values did not make a signifi- 
cant change on the ship location after 600 seconds simula- 
tion, so it can be said that the accuracy of the weighting 
factor is not very important. 

A few simulation runs were performed by changing the 
optimal compensator parameters a small amount and the cost 
curve was plotted. Figure ( 6.10) shows the cost curves for 
three different K values versus Z1 and PI when the ship is 
in calm water, ship speed 15 knots and 1 degree course 
change, and Figure ( 6.11 ) shows the cost curves for ship 
speed 15 knots, encounter angle 60.0 degree, encounter 
frequency 2.5 radian/ second , and sea state 9. 
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lT> 

O 

O 




The optimal parameter values to cost value 
K=0.53, Zl=51.46 , Pl=18.08 and J=0. 0186974 



Figure 6.10 



The Cost Curves vs. Z1 and Pi when Parameters are 
Changing Around the Optimal Values 
for the Calm Water Case. 
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Figure 6.11 The Cost Curves vs. Zl and Pi when Parameters 
Changing Around the Optimal Values 
for the Sea State Case. 
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The results showed that the cost curve is flat around 
the minimum. The same conclusion is made in [Ref. 7 and 8] 
for the SL/7. Using controller parameters on the flat 
portion of the cost curve, but not at the optimum did not 
make a significant differences in simulation results, after 
600 seconds of cruising there was little change in the final 
location of the ship. This property of the cost surface may 
be useful in reducing the time required to find a pratical 
minimum for the cost function. 

The ship speed effects on compensator parameters for 
regular seas can be seen by comparing Tables X, XVI and XVII. 
Observe that for both cases increasing the ship speed made 
the parameter values and the cost value decrease. The same 
observation can be made by comparing the results that are 
found for the SL/7 container ship for different speeds in 
[Ref. 9]. These results strongly indicate that the dynamics 
of the ship determines the optimum structure for the 
controller . 

To see the stability of the system, the ship's third 
order Nomoto model was cascaded with the compensator and the 
open loop system BODE plot was drawn, in every case the 
system is stable. The phase margin varies between 40 
degrees and 70 degrees and the zero cross over of the magni- 
tude curve is around 0.04 radian/ second . It was observed 
that small changes in the compensator parameters do not 
affect stability. With in large limits of parameter varia- 
tion the system is always stable. For regular seas, ship 
speed 13 knots, encounter angle 30.0 degrees, encounter 
frequency 0.6 radian/ second and sea state 6, using compen- 
sator 'A' and 'B', the structure of the system is presented 
in Figure (6.12 ) and the system open loop BODE plot and 
NICHOLS plot are shown in Figures ( 6.13 ) and ( 6.14 ). 
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COMP A 




Where : 





Comp 


. 'A' 


Comp . ' B ' 


Nomoto 


K-- 


-- 0. 


34 


0.30 


PK=0 . 189 


Zl- 


-- 75 


.56 


61.97 


PZ1=18 . 34 


Pl- 


-- 42 


.88 


9.85 


PP1=7 . 67 


P2- 






9.72 


PP2=116 . 93 



Figure 6.12 Open Loop Steering Model. 
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NICHOLS PLOT - GRIN VS PHASE 
o USING ’fl’ TYPE CflSCRDE COMPENSATOR. 




-199.00 -129.00 -120.09 -112.09 -I9<«.09 -99.00 -09.00 -09.00 



Figure 6.13 Open Loop System BODE and NICHOLS Plots. 
^ ^ (Using Comp. ’A' ) . 
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NICHOLS PLOT - GAIN VS PHASE 
o USING ’B’ TYPE CASCADE COMPENSATOR. 




Figure 6.14 Open Loop System BODE and NICHOLS Plots. 

(Using Comp. 'B' ) . 
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Changing the environmentel conditions changes the 
optimal controller parameters. A few simulations were run 
changing the environmental conditions only but the optimal 
controller parameters were kept unchanged The behavior of 
rudder and heading error were observed. Three cases were 
defined depending on the enviromental conditions. 



Case I : Encounter angle 90 degree, encounter 
0.5 radian/ second , sea state 6 and ship speed 
Optimal controller parameters for Case I are 

K = 0.56 Z1 = 39.71 PI = 11.87 

Case II : Encounter angle 90 degree, encounter 
0.5 radian/ second , sea state 8 and ship speed 
Optimal controller parameters for Case II are 



frequency 
1 5 kno t s . 



frequency 
15 knots . 



K = 0.53 Z1 = 27.45 Pi = 08.04 

Case III : Encounter angle 30 degree , encounter frequency 
0.5 radian/ second , sea state 6 and ship speed 15 knots. 
Optimal controller parameters for Case III are 

K = 0.358 Z1 = 66.60 PI = 24.61 

Figure ( 6 . 15 ) and figure ( 6.16 ) show rudder and 
heading error when the ship is in Case I condition using 
Case I and Case II parameters. 

Figure ( 6.17) and figure ( 6.18 ) show rudder and 
heading error when the ship is in Case II condition using 
Case II and Case I parameters. 

Figure ( 6 . 19 ) and figure ( 6.20 ) show rudder and 
heading error when the ship is in Case III condition using 
Case III and Case I parameters. 

Figure ( 6.21 ) and figure ( 6.22 ) show rudder and 
heading error when the ship is in Case I condition using 
Case I and Case III parameters. 
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0.50 -0.25 0.00 0.25 0.50 0.75 
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Figure 6.15 Rudder Motion in Case I with Case I 
and Case II Parameters. 
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Figure 6.16 Heading Error in Case I with Case I 
and Case II Parameters. 




TIME (SEC.) 



Figure 6.17 Rudder Motion in Case II with Case II 

and Case I Parameters. 
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Figure 6.18 Heading Error in Case II with Case II 

and Case I Parameters . 




Figure 6.19 Rudder Motion in Case III with Case III 

and Case I Parameters. 
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TIME (SEC.) 



Figure 6.20 Heading Error in Case III with Case III 

and Case I Parameters. 




Figure 6.21 Rudder Motion in Case I with Case I 
and Case III Parameters. 
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Figure 6.22 Heading Error in Case I with Case I 
and Case III Parameters. 



As seen from Figures (6.15), (6.16), (6.17) and 
( 6.18 ) using Case I optimal parameters in Case II or Case 
II optimal parameters in Case I did not make a big differ- 
ence in rudder motion and heading error, except in the tran- 
sient response part. 

Figures ( 6.19 ) and (6.20 ) show that to use Case I 
optimal parameters when ship is in Case III is not proper 
because those parameters increase rudder and heading error. 

As seen from Figure ( 6.21 ), using Case III parameters 
in Case I provided better rudder motion after the transient 
response part . Cost values were calculated and it was 
observed that when the ship is in Case I using Case I 
optimum parameters it gave a better cost value than Case III 
parameters did, end of the 600 seconds simulation. But it 
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was also observed that if the final time of simulation is 
extended, the difference between the cost values decreases 
and if we continue to increase the final time, using Case 
III parameters in Case I gives a smaller cost value than 
Case I parameters do . 

These results show that the transient response part is 
important in finding the optimum control parameters. 
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VII. AN APPROACH TO AN ADAPTIVE AUTOPILOT 



As seen from the previous chapter the optimal controller 
parameters change for changes in ship conditions such as 
ship speed, loading etc. and also changes in environmental 
conditions such as sea state, encounter angle encounter 
frequency and depth of water. To maintain optimal steering 
performance automatically in the presence of changing condi- 
tions, it is necessary to design an adaptive system which is 
capable of self adjustment of the controller parameters to 
provide minimum added drag due to steering. 

The steering control system, if desired as an adaptive 
autopilot, would consist of four Subsystems as shown in 
Figure ( 7.1 ) . 

• Subsystem #1 would be a computer which will perform to 
find the optimal control parameters. It should get the 
information about ship steering characteristics from 
system state sensors such as a gyrocompass, rudder 
angle potentiometer and speed log and it should feed 
the control parameter values to the controller. 

• Subsystem #2 would be the controller, which should be 
adjustable. It gets the parameter values from 
subsystem #1 and sends the rudder command to subsystem 
#3 which steers the ship. 

• Subsystem #3 is the plant which includes the system 
state sensors and ship steering devices such as servo 
motors, hydraulic pumps, rudder, steering gear etc. 

• Subsystem #4 is a manual control option for safety 
rules . If manual steering is needed it cuts the 
connection between controller and steering devices and 
gives the control to the helmsman. In case of computer 
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failure it cuts the connection between computer and 
controller and sends to the controller parameter values 
which are chosen by the watch officer. It contains 
look up tables which specify control parameters as 
functions of steering and environmental conditions. 




Figure 7.1 Adaptive Control Scheme. 



Success of the this adaptive autopilot system depends on 
the computer program as well as the accuracy of the system 
sensors. The computer program which was used in this 
research may be used on board, but the present program mini- 
mization subroutine needs a lot of computation time and it 
also needs starting guesses for parameters which are 
different for every condition. If computation time is 
reduced to a reasonable time and starting values are made 
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available in proper limits, the modified function minimiza- 
tion program would be considered as an adequate program for 
on board purposes. The work on reducing the computation 
time is presented in [Ref. 7]. 

In the future, ships could have better measurement of 
navigation than can be provided by conventional equipment on 
board. For example, the U.S Navy is involved in a program 
to build a system which is called NAVSTAR/ GLOBAL POSITION 
SYSTEM (GPS). The system will provide extremely accurate 
three-dimensional position and velocity information to users 
anywhere in the world. And also another system is called 
NAVY REMOTE OCEAN SENSING SYSTEM (NROSS) will be able to 
determine wind velocities over th world's oceans with an 
accuracy sufficient to determine ocean surface waves. Using 
such valid information the watch officer can use look up 
tables and insert them into the computer, so system opera- 
tion will be very close to the minimum cost value and the 
function minimization program will accomplish the fine 
tuning rapidly. Detailed information about GPS and NROSS 
can be found in [Ref. 10, 11, 12, 13]. 
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VIII. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE STUDY 



A. CONCLUSIONS 

The conclusions resulting from this research of the 
Mariner class ship fuel consumption as it relate to steering 
might be listed as follows. 

• A steering control system would minimize propulsion 
losses due to steering and maintain desired heading 
with reasonable heading error in every ship and enviro- 
mental condition. 

• The study shows that the best model to describe the 
dynamics of the ship is the Taylor's series expansion, 
which allows both linear and nonlinear terms in the 
ship's equation of motions. Also the third order ship 
Nomoto model is reasonable to use instead of the ship's 
equation of motions. It involves both the sway and yaw 
equations . 

• It is believed that the cost function which is 
presented in Chapter 6 and is commonly used by many 
researchers is an adequate function for on board use. 
It has variables such as heading error and rudder angle 
which can be easily measured on board, and a weighting 
factor (LAMBDA) which is also easy to calculate but 
depends on ship conditions such as ship speed. Results 
in chapter 7 show that accuracy of the LAMBDA value 
does not make significant changes in the controller. 

• In this study two different types of controller were 
tried which have been called controller 'A' and 'B'. 
The structure of these controllers is shown in Chapter 
7. Controller 'A' was determined to be a best struc- 
ture, because in some cases the adaptive calculations 
for controller 'B' did not converge. 
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• An adaptive controller which minimizes propulsion 
losses due to steering is needed when ship characteris- 
tics and environmental conditions change. 

• For performance in fuel saving, an adaptive controller 
may be better than the existing Universal Gyropilot 
(UGP) . 

• Since, equations of motion of surface ships differ only 
in the numerical values of the coefficients, the simu- 
lation programs used for the Mariner class ship would 
be useable for other type of ships, knowing their hull 
characteristics. The studies made for the SL/7 
containership are examples. [Ref. 6, 9, 7, 8]. 

B. RECOMMENDATIONS FOR FUTURE STUDY 

This research does not cover all ship and enviromentel 
conditions, so to get a better understanding about optimal 
controller parameters, it is recommended that the methods be 
applied to find the controller parameters for an expanded 
range of operating conditions. 

This thesis investigated only course_keeping with 
emphasis on minimizing rudder and yawing activity to reduce 
fuel consumption. If a track following control or an auto- 
matic control for replenishment at sea is desired, a 
different cost function might be needed. The nature of the 
cost function for such applications should be studied. 

Irregular seas case were not considered in this study. 
For future work it is necessary to study the effect of 
irregular seas on the controller parameters and on the cost 
function. It is believed that compering the regular sea 
results with irregular sea results would give better under- 
standing about ship's steering characteristics. 
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The transient response has a big affect on finding the 
optimum parameter values. For future studies it may be 
better to increase final time so that the effect of the 
transient response would not be significant. 

Recent studies on roll stabilization shows that using 
rudder stabilization is successful in reducing roll, 
[Ref. 14, 15]. Adding the roll equation into the ship model 
and determining a proper cost function for minimum roll 
motion a controller could be designed for roll 
stabilization. 
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APPENDIX A 



THESIS FORTRAN 



$JOB 



REAL*8 L,L2 .L3 ,L4,L5 ,L6 

REAL*8 X ! XD6t , Y , YDOT , U . UDOT , V , VDOT , YAW , R , RDOT 
REAL-8 TiME,EtlME.Xl,X2,X3,x4,X5 ,X6 ,X7 ,X8 
REAL-8 Y0,Yi,Y2 ,Y^ ,y4,Y^,Y^ ,Y^,Y^ 

-n r' A T -VO XT n ’ M ’ XT O XT O XT /. XT c: XT ^ XT 7 XT Q 



REAL-8 N0,N1,N2,N3 ,N4,N5 ,N6 ,N7,N8 
REAL*8 Cl C2 ! C3 , C4 , C5 , F 1 , F2 , F3 
REAL-'^8 RO,DELT,S,DU,ui,K,Zl,Z2,Pl,P2 
REAL*8 DYAWE , YAWE , YAWC , ISR , ISE , TDIFF , LAMDA 
REAL-8 S1,S2,DS1,DS2,D 
REAL-8 MASS ,IZ ,XG, YVDOT.NVDOT 
REAL*8 YR , YRDOT , NR , NRDOT , FX , FY , MZ 
REAL*8 Rx!rY,RZ,TX,TY,TZ,WA,WE 
REAL - 8 RXR , RYR , RXI , R Y I , MZ R , MZ I 
C INITIAL CONDITIONS FOR INTEGRATION 
C SIMULATION END TIME IN SECONDS 
ETIME=210 . 

TIME=0.0 
I COUNT =1 

C INITIALIZE THE COST FUNCTION 
ISE=0 . 0 
ISR=0.0 
TDIFF=0.0 
LAMDA=2. 91 

C GAIN COEFFICIENTS TO BE OPTIMIZED 
K=0.5 
Zl=43 . 32 
Pl=9 .26 



C X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 

x=o!o 

Y=0 . 0 
XDOT=0 .0 
YDOT=0 .0 



C U, UDOT. V, VDOT ARE FIX COORDINATES ON SHIP 
F1=6.6 
F2=0 . 0 
F3=0.0 
V = 0.0 
UDOT=0.0 
VDOT=0 . 0 
YAW=0.0 
YAWE=0 .0 
R=0 . 0 
RDOT=0 .0 

C ORDERED SPEED IN FEET /SEC 

C 15. *1.689 FT/SEC=15 KNOTS 
Ul=15 .*1. 689 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED 
U=U1 

C D = RUDDER ANGLE 
D=0. 0/57. 296 
L=520. 

L2=L**2 

L3=L*L*L 

L4=L*L3 

L5=L*L4 

L6=L*L5 

C FORCES IN X.Y DIRECTION COMPUTED IN FORCES 

C MOMENTS IN Z 



(UC) 
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5 

3 

5 



[EIGHT: (SEA STATE 2) 



FX=0 . 

FY=0 . 

MZ = 0. 

RXR=- . 37126D5 
RXI= . 68406D5 
RYR=- . 39983D6 
RYI= . 2447D6 
MZR= .296D8 
MZI=- . 17546D8 
RX= (RXR**2+RXI"“2 
RY= (RYR“ “Z + RYI-'-'Z 
RZ = 1 MZR— 2 + MZ I **2 

tx=i)atan2Trxi,rxr; 

TY=DATAN2(RYI ,RYI 
TZ=DATAN2(MZI ,MZR' 

C SIGNIFICANT WAVE 

WA=25 . 

C ENCOUNTER FREQUENCY : (WHEN ENCOUNTER ANGLE IS 00) 

WE=0.75 

C ADDED MASS AND ADDED INERTIA TERMS: 

MASS= . 14685D+07 
IZ= . 23567D+I1 
XG=-23.9 
YR=- . 95066D8 
YRDOT= . I22IID8 
NR=- . 13152DI1 
NRDOT=- . 72177D10 
YVDOT=- . 72459D6 
NVDOT=- . 53846D8 
C 

C HYDRODYNAMIC COEFFICIENTS ARE INSERTED HERE AS PARAMETERS 
RO=I.9876'^5 
YAWE=0.0 
DSI=0.0 
DS2=0.0 
S1=0 .0 
S2 = 0 .0 

200 CONTINUE 
C INPUT YAW COMMAND 
YAWC=0.0 



C 

C 



IF (TIME. GE. 0.0) 
ERROR SIGNAL TO DRI 
YAWE=YAW - YAWC 



YAWC=(1. 0/57. 296) 

VE RUDDER (YAW ACTUAL - YAW ORDERED) 



C 

C 



S=( (;U"U)+ (V-V) ) — .5 
DU=U-U1 

DS1=(YAWE - S1)/P1 
D=(S1 + DSI"Zlj"K 

AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 
Xl= (-0.0253 )*TrO*L2"S) 

X2 = ( 0 . 00 9 48 1 " ( RO ''L2 ) 

X3= (- 0. 00217 )''(R0*L2/S) 

X4=(-0. 189)''(RO*L2) 

X5 =( 0.003 79 }-('RO*l4) 

X6= (-0.02)“ (R0*L2*S"S) 

X7=(0. 168}''UO“L3) 

X8= (0.0196 )“(R0“L2"S) 

LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
Y0=T-0 . 0008 )*(RO*L2“S*S ) 

Y0=0.0 

-0.244)“(RO“L2“S) 



Yl= 

Y2 = 

Y3 = 

Y4 = 

Y4= 

Y5 = 

Y6 = 

Y7 = 

Y8 = 

MOMENT ABOUT 



-1. 702)*(RO*L2/S 
-0.0008)“2r0“L4“o 
-0. 105)“(RO''L3“S) 
YR-MASS) 



^^S) 



3.23)*(RO"L3/S) 

^ r.,L6r.\ '•-7rO"L2 “ S “ S ) 

]"'(RO“L2*S“S) 

Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 



0.0586^ 
-0.00973 
0.23)^'^(R 
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NO = ( 0 . 0 0 0 5 9 ) " ( RO " L 3 “ S •' S ) 
NO=0.O 



Nl = 
N2 = 
N3 = 
N4 = 
N4 = 
N5 = 
N6 = 
N7 = 
N8 = 
COMMON 
Cl = 
C2 = 



-0.0555)“fRO''L3-'S) 
'0. 345)"(RO''L37 s) 

0. 00264 )"(R0-'L3'-' 



-0:0349 j*fR0“L4-S j 



NR-MASS"XC"S) 

-1. 158}*(RO*U,_, 

-0 . 0293)"(R0“L3“8"S j 



0. 00482 )*iR0*L3"S''S 
0. 1032)-(RO*L3) 
COEFFICIENTS 
0.177 
0.327 



C2=(MASS- 



R0*L3 ) 
R0*L3 ) 
OT) 



C3=(-0. 0018)*(RO*L4) 

C3= (mass*xg-yrdotT 

C4=(0.0175}"(RO*L5) 

C4=(IZ-NRD0T) 

C5= No. 00478 j*(R0*L4) 

C 5 = (MA S S "XG - NVDOT ) 
REGULAR WAVES 
FX=WA''RX*DCOS (WE''TIME + TX ' 
FY=WA''RY“DCOS f WE"TIME + TY 
MZ=WA*RZ"DCOS (WE*TIME+TZ 





IF( 


[TIME. EQ 




IF( 


DABS(FY 




IF( 


DABSfMZ 


c 


EQUAT] 


tONS OF 1 




FI 


= XI “DU 




1 


+ X5*R'': 




F2 


= YO + 




1 


+ Y5'“R*' 




F3 


= NO + : 


c 


1 


+ N5"R“’ 



+ X2*DU"DU + X3"DU*DU*DU + X4"V'’'V 
R + X6*D-D + X7-'V'’'R + X8*V->D + FX 
Y1"V + Y2*V*V*V + Y3*V''D*D + Y4-R 
V*Y + Y6-'D + Y7'>D"D"D + Y8''D*V*Y + 
Nl^W + N2“V*V''V + N3'''V*D'>D + N4-'R 
V-V + N6-D + N7-'D"D*D + N8*D"V"V + 



MZ 



UDOT = Fl/Cl 

VDOT = (C4*F2-C3*F3)/ (C2*C4-C5*C3) 
ROOT = fC2*F3-C5*F2)/ N2''C4-C5*C3) 



C WHEN TO PRINTOUT 

IF (TIME. EQ. 0.0) 



GO TO 
GO TO 



IF (ICOUNT 
GO TO 300 
C CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW"57.296 
RDEG = R''57 .296 
RDDEG=RDOT'-57 . 296 
DDEG = D''57 .29 6 
YAWC=YAWC"57.296 
WRITE (6,100) TIME,U,V,R 
100 FORMAT (' ' ,lX,F9.2,li,^’9 

icountN 

C TEST IF WANT TO STOP 

300 IF (TIME.GT.ETIME) GO TO 
C INTEGRATION STEP SIZE DELT 
DELT=1. 

C INTEGRATION 

U=U+UDOT*DELT 

V=V+VDOT“DELT 

R=R+RDOT*DELT 

YAW=YAW+R*DELT 

S1=S1+DS1-DELT 

TIME=TIME+DELT 
ICOUNT=ICOUNT+l 
ISE=ISE + LAMDA"YAWE"*2 
ISR=ISR + D'"‘'2 



50 

50 



6,1X,F9.6,1X,F9.6) 



400 
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TDIFF=ISE+ISR 
GO TO 200 
400 CONTINUE 
STOP 
END 

ENTRY 

END 



APPENDIX B 



DATA FOR SEA STATE PROGRAM 



The sea state program which is explained in [Ref. 4], 
needs information about the ship to calculate the ship's 
added mass, added inertia and seaway disturbance forces and 
moments. Information about the Mariner was gathered from 
[Ref. 27, 28] and presented here in the form which the sea 
state program needs. The current line is drawn on next page 
to show the location of the values in the format 
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APPENDIX C 



PROGRAM TO CALCULATE OPTIMAL GAINS 



Algorithm used here is showed in Figure ( C.l ). 




Figure C.l Algorithm to Calculate Optimal Gains. 



The disturbance forces (FX,FY) and moment (MZ) for 
regular seas were aplied in program as a cosine wave. The 
procudure to do this is: 
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• The real and imaginary value of force or moment is read 
from the sea state program output depending on ship 
speed, encounter frequency and encounter angle. 

• These values are converted to magnitude and phase 
values . 

• Using the formula below, forces and moment are created 
and added into the ship's equations of motion. 

Force or Moment = WA“MAGNITUDE-cos (WE-TIME + PHASE) 

Where : 

WE = Encounter frequency (radian/second) 

WA = Significant wave height (ft) 

The correspondence between sea state and significant wave 
height is indicated in Table XXI [Ref. 29]. During this 
research the values that are presented in Table XXII were 
used as significant wave height (WA) . 



TABLE XXI 

Sea State vs Range for Wave Height 



Sea State 
(Beauf oFE s^ale ) 



Range for .wave . height 
(Feet ) 



0 

1 

2 

3 

4 

5 

6 

7 

8 
9 




3.28-4.92 

6.56-8.20 

9.84-13.1 

13.1-18.2 



18.2-24.6 

23.0-32.9 
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TABLE XXII 

Sea State vs Wave Height 



Sea state 

5 

6 

7 

8 
9 



Wave height 



5 

10 

17 

20 



0 

0 

5 

0 



25.0 



(ft) 



The program is set up to calculate the optimal gains for 
controller A. It can be modified to obtain optimal gains for 
the rest of the controllers. After obtaining the optimal 
gains the program must be modified to do a simulation. The 
program has sufficient comments for appropriate changes. 

This program can be modified to obtain the Nomoto models. 
It is referenced in Chapter 4. The following need to be 
changed . 

C GAIN COEFFICIENTS TO BE OPTIMIZED 
K=XX(1) 

Z1=XX(2) 

P1=XX(3) 

P2=XX(4) 

C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL - YAW COMMAND) 
c FOR EQUATIONS OF MOTION. 

D=YAW - YAWC 

C ERROR SIGNAL TO DRIVE RUDDER (YAW COMMAND - YAW ACTUAL) 

C FOR NOMOTO 3RD ORDER MODEL. 

D2=YAWC-YAW2 
DQ1=(D2 - Q1)/P1 
DQ2= ( (K-’'( (Z^’^DQ1) + Q1)/P2 
C INTEGRATION 

Q1=Q1+DQ1"DELT 

Q2=Q2^DQ2"DELT 
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YAW2 = YAW2 + Q2 •' DELT 
COST FUNCTION 

TDIFF=TDIFF+ (YAW- YAW2 )**2 
PROGRAM TO CALCULATE OPTIMAL GAINS FOR CONTROLLER 



//TANSAN JOB ( 1789 , 0356 RESEARCH CLASS=J 
//“MAIN 0RG=NPGVM1. 1789P 
// EXEC FRTXCLGP,IMSL=DP,REGION=1024K 
//FORT. SYS IN DD “ 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 



C OBTAINED 
XLC^-^) 



CHANGE XS(“) TO X(“) AND DELETE XU("),AND 




C 

C 



DIMENSION XS(3) ,XU(3) ,XL(3) 

XS(1)=0.6 

XSU) = 41.58 

XS(3)=11. 33 

I) IS THE STARTING GUESS 
I) IS THE LOWER LIMIT FOR THE 
I) IS THE UPPER LIMIT FOR THE 




I’TH VARIABLE 
I’TH VARIABLE 



C 

C 

C 



25 

30 

40 



C 

C 



= . 3 
= 0.9 
= 30. 

= 50. 

= 6. 

= 15. 

TION OF THE FOLLOWING PARAMETERS 
IS DISCUSSED IN BOXPLX 
R=9. /13. 

NTA=1000 
NPR=0 
NAV=0 
NV=3 
IP = 0 

THE FOLLOWING STATEMENT MUST BE CHANGED TO 
CALL PLANT (_XX) 

IF ONLY SIIWLATION IS WANTED 

CALL BOXPLX? NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , lER ) 
WRITE (6,25) 

FORMAT (li,’ OPTIMAL GAINS',/) 

DO 30 I=i!3 
WRITE(6 ,40)1 ,XS(I) 

FORMAT UX, 'X? ' ,12 , )= ' ,F14. 7) 

STOP 

END ^ ^ 

SUBROUTINE PLANT (XX) 

SUBROUTINE PLANT ( XX ) SIMULATES THE SHIP 
COMMON TDIFF 
REAL“8 L,L2 ,L3 ,L4,L5 ,L6 

REAL" 8 X , XD6t , t , Yf)OT , U . UDOT , V , VDOT , YAW , R , RDOT 

REAL“8 TIME,EtlME.Xl,X^,X3,X4,X5 ,X6,X7,X8 

REAL*8 Y0,y1,Y2,Y^,y4,y3 ,Y^ ,Y^ ,Y8 

REAL*8 NO , N1 , N2 , N3 , N4 , N5 , N6 , N7 , N8 

REAL'-^8 Cl C2,C3,C4,C5,F1,F2,F3 

REAL*8 RO,DEtT,§,DU,ui,K,Zl,Z2,Pl,P2 

REAL "8 DYAWE,YAw^;,YAwC,ISR,iSE,TDIFF,LAMDA 

REAL“8 SI ,S2 ,DS1,DS2 ,D 

REAL ''8 MA^S , iZ ,XG , YVDOT .NVDOT , YR,YRDOT 

REAL* 8 NR , NRDOT , FX , FY , MZ , RXR , RYR , RXI , RYI , MZR , MZI 

REAL*8 RX , RY , RZ , TX , TY , TZ , WA , WE 

DIMENSION XX(3) 

INITIAL CONDITIONS FOR INTEGRATION 
SIMULATION END TIME IN SECONDS 
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on 



c 

c 

c 

c 



c 

c 

c 

c 



c 

c 

c 



ETIME=600 . 

TIME=0. 0 
IC0UNT=1 

INITIALIZE THE COST FUNCTION 
ISE=0 . 0 
ISR=0.0 
TDIFF=0.0 
LAMDA=2 .91 

GAIN COEFFICIENTS TO BE OPTIMIZED 
K=XX(1} 

Z1=XX(2) 

P1=XX( 3 ) 

X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0.0 
XDOT=0.0 
YDOT=0.0 



U,UDOT.V,VDOT ARE FIX COORDINATES ON SHIP 
F1=6.6 
F2=0.0 
F3=0.0 
V = 0.0 
UD0T=0.0 
VD0T=0.0 
YAW=0.0 
YAWE=0.0 
R=0.0 
RDOT=0.0 

ORDERED SPEED IN FEET/ SEC 
15. “1.68 9 FT/SEC=15 KNOTS 
Ul=15 .^n. 689 

AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=U1 

D = RUDDER ANGLE 
D=0. 0/57 .296 
L=520. 

L2 = L"''2 
L3 = L*L''L 
L4=L“L3 
L5=L*L4 
L6=L*L5 

FORCES IN X.Y DIRECTION COMPUTED IN FORCES 

MOMENTS IN Z 

FX=0. 

FY=0, 

MZ=0. 

RXR= . 33374D3 

RXI= . 16064D3 

RYR=- . 17254D5 

RYI= . 1435D4 

MZR= .30976D7 

MZI= .84672D6 

RX= (RXR’’^*2 + RXI '"’'2 )*“ . 5 

RY= ( RYR “ *2 + RYI **2 ) “ “ . 5 

RZ= tMZR*“2 + MZI*"2 P— . 5 

TX=DATAN2(RXI,RXR) 

TY=DATAN2 ( RYI , RYI ) 

T7=DATAM? ( M7T M7P i 

SIGNIFICANT wAve Height state 2) 

WA=25. 

ENCOUNTER FREQUENCY : (WHEN ENCOUNTER ANGLE IS 00) 
WE=2.5 V J 

ADDED MASS AND ADDED INERTIA TERMS: 

MASS= . 14685D+07 
IZ= .23567D+11 
XG=-23 . 9 
YR=- .24965D8 
YRDOT= . 75199D7 
NR=- .47281D10 
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NRDOT=- . 29233D11 
YVDOT=- . 19723D6 
NVDOT=- .23146D8 



INSERTED HERE AS PARAMETERS 



HYDRODYNAMIC COEFFICIENTS ARE 
R0=1. 9876^^5 
YAWE=0. 0 
DS1=0. 0 
DS2=0.0 
S1=0.0 
S2=0.0 

200 CONTINUE 
INPUT YAW COMMAND 
YAWC=0.0 

IF (TIME. GE. 0.0) YAWC= ( 1 . 0/ 57 . 296 ) 

ERROR Signal to drive rudder(yaw actual 

YAWE=YAW - YAWC 
S = ( 7u "U ) + ( V * V ) ) ** . 5 
DU=U-U1 

DS1=(YAWE - S1)/P1 
D=(si + DSl'-ZlI-'K 

AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 
Xl= ( - 0 . 025 3 ) -''7 rO*L2“S ) 

X2=(0. 00948 j*Uo-L2) 

X3= ( - 0 . 002 17 ) (R0"L2 / S ) 

X4=(-0. 189)”(rO^‘'L2) 

X5= (0.0037 9 )''(R0*l4) 

X6 = ( - 0 . 02 ) * ( R0'H.2*S - S ) 

X7=(0. 168l''(RO*L3) 

X8=(0.0196)" (RO *L2 -'S ) 

LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
Y0=(-0. 0008)*(R0-L2*S*S) 

Y0 = 0.0 



YAW ORDERED) 



Yl = (-0. 244)" (R0''L2-S) 
Y2=(-l. 702j*(RO"L2/S) 
Y3 = ( - 0 . 0 0 0 8 ) " 7 R0'-’^L4* S 
Y 4 = ( - 0 . 1 0 5 ] " ( R 0 * L 3 " S ) 
YR-MASS) 

*(RO" 



^S) 



Y4 = 

Y5 = 

Y6 = 

Y7 = 

Y8 = 

MOMENT ABOUT 



3V23)*(RO*L3/S) 

' 0586 )"(RQ*L2"S"S2 



-0. 00975]'MRO*L2*S 
0.25)*^Rp"L2) 



S) 



AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 



N1 
N2 = 
N3 = 
N4 = 
N4 = 
N5 = 
N6 = 
N7 = 
N8 = 
COMMON 
Cl = 



N 0 = ( 0 . 0 0 0 5 9 ) " ( RO * L 3 S " S ) 
N0 = 0 ^ 



-0.0555)"(RO"L3"S) 

'0.345)"(RO"L37s) 

' 0 . 0 0 2 6 4 ) ■" 7 RO*L 3 " S ) 
■-0.0349)"'(RO*L4"S) 
'NR-MASS"XG"S) 

-1. 158)"(RO"L4/S) 

' - 0 . 0 2 9 3 ) " ( RO " L 3 S S ) 
0.00482)*] RO*L 3 * S * S ) 
0. 1032)"(RO*L3) 



COEFFICIENTS 
__ ,0. 177)*(RO*L3) 
C2=(0.327j*]RO*L3) 



C2= (MASS-tVOOT) 
C3=(-0.0018 )*]rO*L4) 

C3= (MASS'-'XG-YRDOT) 
C4=(0.0175)*(RO*L5) 

C4= (IZ-NRDOt) 

C5=(-0. 00478 j*(R0*L4) 
C5=(MASS*XG-NVDOT) 
REGULAR WAVES 
FX=WA*RX*DCOS ( WE'"TIME+TX , 
FY=WA*RY*DCOS ( WE*TIME+TY , 
MZ=WA*RZ*DCOS ( WE-'TIME + TZ , 
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IF(TIME.EQ.O.O) FX=0.0 ^ ^ 

IF (DABS (FY) .LT. 0.00000001) FY=0.0 
IF (DABS (MZ) .LT. 0.00000001) MZ=0.0 
EQUATIONS OF MOTION . . , . , . , 

FI = X1*DU + X2-'DU*DU + X3*DU*DU“DU + X4“V"V 
+ X5“R''R + X6*D*D + X7*V"R t XS^V-'D + FX 
F2 = YO + Y1"V + Y2*V-V''V + Y3*V"D"D + Y4"R 
+ Y5-’'R*V*V + Y6“D + Y7*D"D*D + Y8*D“V'’'Y + 
F3 = NO + Nl-W + N2'>V''V"V + N3*V*D*D + N4"R 
+ N5“R"V"V + N6*D + N7'''D''D''D + N8“D"V*V + 



FY 

MZ 



= (c4''F2-C3*F3)/ (C2“C4-C5*C3 
= (C2*F3-C5“F2)/ (C2“C4-C5*C3 



UDOT = Fl/Cl 
VDOT - 
RDOT 

C WHEN TO PRINTOUT 

IF TiCOUNT. EQ. 21) GO TO 50 
GO TO 300 

C CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW''57.296 
RDEG=R*57 .296 
RDDEG=RDOT*57 .296 
DDEG=D"57 .296 
YAWC=YAWC-57 . 296 
WRITE (6,100) TIME.TDIFF 
100 FORMAT (’ ’ ,1X,F10.2,1X,F20.10) 

ICOUNT=l 

C TEST IF WANT TO STOP 

■ 300 IF (TIME.GT.ETIME) GO TO 400 

C INTEGRATION STEP SIZE DELT 
DELT=1. 

C INTEGRATION 

U=U+UDOT"DELT 

V=V+VDOT'-DELT 

R=R+RDOT“DELT 

YAW=YAW+R''DELT 

S1=S1+DS1*DELT 

S2=S2+DS2*DELT 



C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



TIME=TIME+DELT 
ICOUNT=ICOUNT+l 
ISE = ISE + LAMDA"YAWE — 2 
ISR=ISR + D**2 
GO TO 200 
400 TDIFF=ISE+ISR 

WRITE(6,500) TIME.TDIFF 
500 FORMAT ( ' ' , IX , F9 . 2 , IX , F20 . 10 ) 

RETURN 

END 

DELETE ALL THE FOLLOWING SUBROUTINE IF SIMULATION ONLY 
AND NOT OPTIMIZATION IS WANTED 



SUBROUTINE BOXPLX 
PURPOSE 



(CATEGORY HO) 



BOXPLX IS A SUBROUTINE USED TO SOLVE THE PROBLEM OF 
LOCATING A MINIMUM (OR MAXIMUM) OF AN ARBITRARY OBJECT- 
IVE FUNCTION SUBJECT TO ARBITRARY EXPLICIT AND/OR 
IMPLICIT CONSTRAINTS BY THE COMPLEX METHOD OF M.J. BOX. 
EXPLICIT CONSTRAINTS ARE DEFINED AS UPPER AND LOWER 
BOUNDS ON THE INDEPENDENT VARIABLES IMPLICIT CONSTRAINTS 
MAY BE ARBITRARY FUNCTION OF THE VARIABLES. TWO FUN- 
CTION SUBPROGRAM TO EVALUATE THE OBJECTIVE FUNCTION AND 
IMPLICIT CONSTRAINTS, RESPECTIVELY. MUST BE SUPPLIED 
BY THE USER 7 SEE EXAMPLE BELOWj . b6xPLX ALSO HAS tHE 
OPTION TO PERFORM INTEGER PROGtoOMING, WHERE THE VALUES 
OF THE INDEPENDENT VARIABLES ARE RESTRICTED TO INTEGERS. 
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C USAGE 

C 

C CALL BOXPLX (NV , NAV , NPR , NTA, R ,XS , IP , XU , XL , YMN , lER) 

C 

c 

C DESCRIPTION OF PARAMETERS 

C 

C NV AN INTEGER INPUT DEFINING THE NUMBER OF INDEPENDENT 
C VARIABLES OF THE OBJECTIVE FUNCTION TO BE MINIMIZED. 

C NOTE: MAXIMUM NV + NAV IS PRESENTLY 50. MAXIMIM NV IS 

C 25. IF THESE LIMITS MUST BE EXCEEDED, PUNCH A SOURCE 

C DECK IN THE USUAL MANNER, AND CHANGE THE DIMENSION 

C STATEMENTS. 

C 

C NAV AN INTEGER INPUT DEFINING THE NUMBER OF AUXILIARY var 
C iaBLES THE USER WISHES TO DEFINE FOR HIS OWN CONVENIENCE. 
C TYPICALLY HE MAY WISH TO DEFINE THE VALUE OF EACH IMPLICI 
C CONSTRAINT FUNCTION AS AN AUXILIARY VARIABLE. IF THIS 
C IS DONE, THE OPTIONAL OUTPUT FEATURE OF BOXPLX CAN BE 
C USED TO OBSERVE THE VALUES OF THOSE CONSTRAINTS AS THE 
C SOLUTION PROGRESSES. AUXILIARY VARIABLES, IF USED, 

C SHOULD BE EVALUATED IN FUNCTION KE (DEFINED BELOW). 

C NAV MAY BE ZERO. 

C 

C NPR INPUT INTEGER CONTROLLING THE FREQUENCY OF OUTPUT 
c desired for diagnostic purposes. 

C IF NPR .LE. 0, NO OUTPUT WILL BE 

C PRODUCED BY BOXPLX. OTHERWISE, THE CURRENT COMPLEX OF 
C K= 2"NV VERTICES AND THEIR CENTROID WILL BE OUTPUT AFTER 
C EACH NPR PERMISSIBLE TRIALS. THE NUMBER OF TOTAL TRIALS, 
C NUMBER OF FEASIBLE TRIALS, NUMBER OF FUNCTION EVALUATIONS 
C AND NUMBER OF IMPLICIT CONSTRAINT EVALUATIONS ARE IN- 
C CLUDED IN THE OUTPUT. 

C ADDITIONALLY, (WHEN NPR .GT. 0) THE SAME INFORMATION 
C WILL BE OUTPUT: 

C 

C 1) IF THE INITIAL POINT IS NOT FEASIBLE, 

C 2) AFTER THE FIRST COMPLETE COMPLEX IS GENERATED, 

C 3) IF A FEASIBLE VERTEX CANNOT BE FOUND AT SOME TRIAL, 

C 4) IF THE OBJECTIVE VALUE OF A VERTEX CANNOT BE MADE 
C NO-LONGER-WORST. 

C 5) IF THE LIMIT ON TRIALS (NTA) IS REACHED AND, 

C 6) WHEN THE OBJECTIVE FUNCTION HAS BEEN UNCHANGED FOR 
c 2"NV TRIALS, INDICATING A LOCAL MINIMUM HAS BEEN 
C FOUND. 

C 

C IF THE USER WISHES TO TRACE THE PROGRESS OF A SOLUTION, 

C A CHOICE OF NPR = 25, 50 OR 100 IS RECOMMENDED. 

C 

C NTA INTEGER INPUT OF LIMIT ON THE NUMBER OF TRIALS 
c allowed in the calculation. 

C IF THE USER INPUTS NTA .LE. 0, A default 
C VALUE OF 2000 IS USED. WHEN THIS LIMIT IS REACHED 
c CONTROL RETURNS TO THE CALLING PROGRAM WITH THE BEST 
C ATTAINED OBJECTIVE FUNCTION VALUE IN YMN, AND THE BEST 
C ATTAINED SOLUTION POINT IN XS . 

C 

C R A REAL NUMBER INPUT TO DEFINE THE FIRST RANDOM NUMBER 
C USED IN DEVELOPING THE INITIAL COMPLEX OF 2"NV VERTICIES. 
C (0. .GT. R .LT. 1.) IF R IS NOT WITHIN THESE BOUNDS, 

C IT WILL BE REPLACED BY I./3. . 

C 

C XS INPUT REAL ARRAY DIMENSIONED AT LEAST NV+NAV. 
c the first nv must contain a 
C FEASIBLE ORIGIN FOR STARTING THE CAL- 

C CULATION. THE LAST NAV NEED NOT BE INITIALIZED. UPON 
C RETURN FROM BOXPLX, THE FIRST NV ELEMENTS OF THE ARRAY 
C CONTAIN THE COORDINATES OF THE MINIMUM OBJECTIVE 
C function, AND THE REMAINING NAV (NAV . GE . 0) CONTAIN THe 
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C values of THE CORRESPONDING AUXILIARY VARIABLES. 

C 

C IP INTEGER INPUT FOR OPTIONAL INTEGER PROGRAMMING. 

C if ip=l, THE VALUES OF THE INDEPENDENT VARIABLES WILL 
C be replaced WITH INTEGER VALUES (STILL STORED AS REAL''4 ) 
C 

C XU A REAL ARRAY DIMENSIONED AT LEAST NV INPUTTING THE 
C upper BOUND ON EACH INDEPENDENT VARIABLE, (EACH EXPLICIT 
C conSTRAINT). INPUT VALUES ARE SLIGHTLY ALTERED BY BOXPLX 
C 

C XL A REAL ARRAY DIMENSIONED AT LEAST NV INPUTTING THE 
c lower bound on each independent 
C VARIABLE, (EACH EXPLICIT CONstraint). 

C NOTE: FOR BOTH XU AND XL CHOOSE REASONABLE 

C VALUES IF NONE ARE GIVEN, NOT VALUES WHICH ARE 
C magnitudes ABOVE OR BELOW THE EXPECTED SOLUTION. 

C input values are SLIGHTLY ALTERED BY BOXPLX. 

C 

C YMN THIS OUTPUT IS THE VALUE (REAL*4) OF THE OBJECTIVE 
C funcTION, CORRESPONDING TO THE SOLUTION POINT OUTPUT IN X 
C 

C lER INTEGER ERROR RETURN. TO BE INTERROGATED UPON 
C return FROM BOXPLX. lER WILL BE ONE OF THE FOLLOWING: 

C 

C =-l CANNOT FIND FEASIBLE VERTEX OR FEASIBLE CENTROID 
C AT THE START OR A RESTART (SEE 'METHOD' BELOW). 

C =0 FUNCTION VALUE UNCHANGED FOR 'N' TRIALS. (WHERE 
C N=6'’-NV+10} THIS IS THE NORMAL RETURN PARAMETER. 

C =1 CANNOT DEVELOP FEASIBLE VERTEX. 

C =2 CANNOT DEVELOP A NO -LONGER- WORST VERTEX. 

C =3 LIMIT ON TRIALS REACHED. (NTA EXCEEDED) 

C NOTE: VALID RESULTS MAY BE RETURNED IN ANY OF THE 

C ABOVE CASES. 

C 

C EXAMPLE OF USAGE 
C 

C THIS EXAMPLE MINIMIZES THE OBJECTIVE FUNCTION SHOWN IN 
C the EXTERNAL FUNCTION FE(X). THERE ARE TWO INDEPENDENT 
C varlABLES X(l) & X(2), AND TWO IMPLICIT CONSTRAINT 
C function X(3) & X(4) miCU ARE EVALUATED AS AUXILIARY 
C variables (see EXTERNAL FUNCTION KE(X) ). 

C DIMENSION XS(4) ,XU(2) ,XL(2) 



CC 

CC 



STARTING GUESS 



C 

C 

CC UP 

c 

CC LO 
C AJL 

C XL 

CC 




AL(1) = 0. 

XL(2) = 0. 

R = 9./ 13. 
NTA = 5000 



IMITS 
= 0.0 
= 0.0 



C R 

C N' 



C 

c 



NPR = 50 



C NV 

C IP 



NAV = 2 
NV = 2 
IP = 0 



CC 

C 

C 

C 



WRIT 

IFORM 

2,A(e 



CALL 




NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , lER ) 
ll,i=1.4),YMN IEA) 

THE POINT IS LOCATED AT (XS(I)=) ' 

TION VALUE IS ’,E13.7,' lER = ’,15) 



C 

CC 

C 

C 



STOP 

END 
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cc 

cc 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

cc 

c 

c 

c 

cc 

cc 

c 

c 

cc 

cc 

c 

c 

c 

€ 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



FUNCTION KEfX) 

EVALUATE CONSTRAINTS. SET KE=0 IF NO IMPLICIT CONSTRAINT 
is VIOLATED, OR SET KE=1 IF ANY IMPLICIT 
constraint is violated. 

DIMENSION X(4) 

XI = X(lj 



.GT. 6 . ) GO TO 1 



.on 



X2 = X(2' 

KE = 0 
x 73) = XI + 1.732051*X2 
IF, (X 3) .LT. 0. .OR. X(3) 
Xl4) = Xl/1. 732051 -X2 
IF (X(4) .GE. 0.) RETURN 

KE = 1 
RETURN 
END 



FUNCTION FE(X) 

DIMENSION X(4) 

THIS IS THE OBJECTIVE FUNCTION. 

FE= -(X(2)“*3 *(9.-(X(l)-3. )**2)/(46. 76538)) 

RETURN 
END 

METHOD 

THE COMPLEX METHOD IS AN EXTENSION AND ADAPTION OF 
the simple method of linear programming. 

STARTING WITH ANY ONE feasible point in n-dimensic 
A COMPLEX OF 2''N vertices is constructed by 
SELECTING RANDOM POINTS WITHIN THE feasible 
REGION. FOR THIS PURPOSE N COORDINATES ARE FIRST 
RANDOMLY CHOSEN WITHIN THE SPACE BOUNDED BY EXPLICIT CON- 
STRAINTS. THIS DEFINES A TRIAL INITIAL VERTEX, 
it is then checked for possible violation 
OF IMPLICIT constraints: if one or more are violated, 

THE TRIAL INITIAL VERTEX IS DISPLACED half of its 
DISTANCE FROM THE CENTROID OF PREVIOUSLY SELECTED initial 
VERTICES. IF NECESSARY THIS DISPLACEMENT PROCESS IS 
REPEATED UNTIL THE VERTEX HAS BECOME FEASIBLE. IF THIS 
fails to happen after 5"'n+10 displacements, 

THE SOLUTION IS ABANDoned . after each vertex is added 
TO THE COMPLEX, THE CURRENT centroid is checked for 
FEASIBILITY. IF IT IS INFEASIBLE, the last trail 
VERTEX IS ABANDONED AND AN EFFORT TO GENERATE an alter- 
ATIVE TRIAL VERTEX IS MADE. IF 5-N+10 VERTICES ARE 
ABANDONED CONSECUTIVELY, THE SOLUTION IS TERMINATED. 

IF AN INITIAL COMPLEX IS ESTABLISHED, THE BASIC 
computation loop is initiated. 

THESE INSTRUCTIONS FIND THE CURRENT WORST vertex, that 
IS, THE VERTEX WITH THE LARGEST CORRESPONDING value for 
THE OBJECTIVE FUNCTION, AND REPLACE THAT VERTEX BY 
ITS OVER- REFLECTION THROUGH THE CENTROID OF ALL OTHER 
vertices, (if the vertex to be 

REPLACED IS CONSIDERED AS A VECTOR IN n-space, 

ITS OVER-REFLECTION IS OPPOSITE IN DIRECTION, IN- 
CREASED IN LENGTH BY THE FACTOR 1.3, AND COLLINEAR WITH 
THE REPLACED VERTEX AND CENTROID OF ALL OTHER VERTICES.) 

WHEN AN OVER- REFLECTION IS NOT FEASIBLE OR REMAINS 

worst, it is considered not-permissible 

AND IS DISPLACED HALFWAY TOWARD the centroid. 

AFTER FOUR SUCH ATTEMPTS ARE MADE UNSUCCESSFULLY 
EVERY FIFTH ATTEMPT IS MADE BY REFLECTING THE OFFENDING 
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c vertex through the present best 

C VERTEX, INSTEAD OF THROUGH THE CENtroid. if 5“'n+10 
C DISPLACEMENTS AND OVER- REFLECTIONS OCCUR without a 
C SUCCESSFUL (PERMISSIBLE) RESULT, THE CURRENT BEST 
C VERTEX IS TAKEN AS AN INITIAL FEASIBLE POINT FOR A 
c restart run of the complete process. . 

C RESTARTING IS ALSO UNDERTAKEN when 6"'nv+10 consecutive 
C TRIALS HAVE BEEN MADE WITH NO SIGNIFicant change in the 
C VALUE OF THE OBJECTIVE FUNCTION. IN ALL cases, 

C RESTARTING IS INHIBITED IF THE LAST RESTART DID NOT 
C PRODUCE A SIGNIFICANT IMPROVEMENT IN THE MINIMUM 
c attained. 

C 

C IT IS RECOMMENDED THAT THE USER READ THE REFERENCE FOR 
C FURTHER USEFUL INFORMATION. IT SHOULD BE NOTED THAT THE 
C ALGORITHM DEFINED THERE HAS BEEN ALTERED TO FIND THE 
C CONSTRAINED MINIMUM, RATHER THAN THE MAXIMUM. 

C 

C 

C 

C REMARKS 

C 

C THE INTEGER PROGRAMMING OPTION WAS ADDED TO THIS PROGRAM 
C AS SUGGESTED IN REFERENCE (2). A MIXED 
c integer/ cont inuous variable version of boxplx 
C WOULD BE EASY TO CREATE BY DEclaring "ip" to be an array 
C OF NV CONTROL VARIABLES WHERE IP (il=l would indicate 
C THAT THE I-TH VARIABLE IS TO BE CONFINED to integer 
C VALUES. EACH STATEMENT OF THE FORM 'IF (IP .EQ.i)' etc. 

C WOULD THEN NEED TO BE ALTERED TO 'IF (IP (I) -EQ. 1)' etc. 
C , WHERE THE SUBSCRIPT IS APPROPRIATELY CHOSEN. NORMALLY, 
C XU AND XL VALUES ARE ALTERED TO BE AN EPSILON 'WITHIN' 
c actual values 

C DECLARED BY THE USER. THIS ADJUSTMENT IS NOT MADE 
C WHEN IP=1. 

C 

C NOTE: NO NON-LINEAR PROGRAMMING ALGORITHM CAN GUARANTEE 

c that the answer found is the global 

C MINIMUM, RATHER THAN JUST A local minimum, however, 

C ACCORDING TO REF. 2, THE COMPLEX method has an advantage 
C IN THAT IT TENDS TO FIND THE GLOBAL minimum more 
C FREQUENTLY THAN MANY OTHER NON-LINEAR PROGRAM- 
C MING ALGORITHMS. 

C 

C IT SHOULD BE NOTED THAT THE AUXILIARY VARIABLE FEATURE 
c can also be used to deal with 

C PROBLEMS CONTAINING EQUALITY CONstraints. any equality 
C CONSTRAINT IMPLIES THAT A GIVEN VARiable is not truly 
C INDEPENDENT. THEREFORE, IN GENERAL, ONE variable 
C INVOLVED IN AN EQUALITY CONSTRAINT CAN BE RENUMBERED from 
C THE SET OF NV INDEPENDENT VARIABLES AND ADDED TO THE SET 
C OF NAV AUXILIARY VARIABLES. THIS USUALLY INVOLVES 
C renumbering THE INDEPENDENT VARIABLES OF THE GIVEN 
C problem 

C SUBROUTINES AND FUNCTIONS REQUIRED 

C SUBROUTINE 'BOUT' AND FUNCTION ' FBV ' ARE INTEGRAL 
C parts of THE BOXPLX PACKAGE. 

C 

C TWO FUNCTIONS MUST BE SUPPLIED BY THE USER. THE FIRST, 
c ke(x), is used to evaluate the implicit 

C CONST^INTS. set KE=0 at the beginning of the function 
C , THEN EVALUATE THE IMPLICIT CONstraints. in the example 
C ABOVE, THE FIRST CONSTRAINT, X(3J, must be within the 
C RANGE (0. . LE . X(3) . LE . 6.J. TH^l SECOND constraint x(4) 
C , MUST BE .GE. 0. . IF EITHER CONSTRAINT IS not within 
C THESE BOUNDS, CONTROL IS TRANSFERRED TO STATEMENT 1, 

C AND KE IS SET TO "1" AND CONTROL IS RETURNED TO BOX^>LX. 

C 
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THE SECOND FUNCTION THE USER MUST PROVIDE EVALUATES THE 
objective function, it is 

CALLED FE(X) AS SHOWN IN THE EXAMple above, and fe 
MUST BE SET TO THE VALUE OF THE OBJECTIVE function 
CORRESPONDING TO CURRENT VALUES OF THE NV INDEPENDENT 
VARIABLES IN ARRAY 'X' . 
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SUBROUTINE BOXPLX (NV ,NAV , NPR , NTZ ,RZ ,XS , IP , BU , BL , YMN , lER) 

DIMENSION V(50 ,50) , FUN(50), SUM(25), CEN(25), XS(NV) 
1 , bu (nv ) , bl (NV ) 

KV = 5 
EP = l.E-6 
NTA = 2000 

IF (NTZ.GT.O) NTA = NTZ 
R = RZ 

IF (R.LE.O. .OR.R.GE. 1. ) R=l./3. 

NVT = NV+NAV 

TOTAL VARS, EXPLICIT PLUS IMPLICIT 

NT = 0 

CURRENT TRIAL NO. 

NPT = 0 

CURRENT NO. OF PERMISSIBLE TRIALS 
NTFS = 0 

CURRENT NO. OF TIMES F HAS BEEN ALMOST UNCHANGED 
CHECK FEASIBILITY OF START POINT 



DO 4 1=1, NV 
VT = xs(i) 

IF (BL7i) .LE. VT) GO TO 1 
II = -I 
VT = BL(I) 

GO TO 2 

IF (BU(I) .GE.VT) GO TO 3 
II = I 



2 

3 



VT = 


BU ( I ) 




(6,49) 




IF (NPR.GT.O) 


WRITE 


II 


v(i ■ 


L) = VT 






CEN(] 


[j = VT^ 








IF (j 


[P.EQ.l) 


GO TO 4 




BL(I 


! = BLti) 


+AMAX1I 


EP,EP* 


ABS 


BU(I' 


) = BU(I) 


-AMAXll 


(ep,ep* 


ABS 


SUM(: 


[) = VT 












89 



nn n nn n rin nn n on on on 



C 

C 



NCE = 1 

NUMBER OF CONSTRAINT EVALUATIONS 

1 = 1 

IF (KE(V(1,1)) -EQ-0) GO TO 5 
IF INPR.LE.O) GO TO 12 
WRITE (6,50) 

GO TO 12 
NFE = 1 



OF VARIABLES. 



NUMBER OF VERTICES (K) = 2 TIMES NO. 

K = 2*NV 

NUMBER OF DISPLACEMENTS ALLOWED. 

NLIM = 5*NV+10 

NUMBER OF CONSECUTIVE TRIALS WITH UNCHANGED FE TO 
terminate . 

NCT = NLIM+NV 
ALPHA = 1.3 
FK = K 
FKM = FK-1. 

BETA = ALPHA+1. 

INSURE SEED OF RANDOM NUMBER GENERATOR IS ODD. 

IQR = R--1.E7 

IF (M0D(IQR,2) .EQ.O) IQR=IQR+101 

SET UP INITIAL VERTICES 
FUN(l) = FE(V(1,D) 

YMN = FUN(l) 

6 FI = 1. 

FUNOLD = FUN(l) 

DO 15 1=2, K 
FI = FI+1. 

LIMT = 0 

7 LIMT = LIMT+1 

END CALCULATION IF FEASIBLE CENTROID CANNOT BE FOUND. 
IF (LIMT. GE. NLIM) GO TO 11 

DO 8 J=1,NV 

RANDOM NUMBER GENERATOR (RANDU) 

IQR = IQR--655 39 

IF (IQR.LT.O) IQR = IQR+2147483647+ 1 
RQX = IQR 

RQX = RQX'’^4656613E-9 

VCJ.I) = BL(Jj+RQX"(BU(J)-BL(J) ) 

IF llP.EQ.l) VP,I)=AINt(v(J,I)+.5) 

8 CQNTINUE 

DQ 10 L=1,NLIM 
NCE = NCE+1 

IF (KE(V(1,I) ) .EQ.O) GO TO 13 



DO 9 J=1,NV 
VT = . 5'-(V( J ,I) + CEN(J 
IF (IP.EQ ” 

V(J,l) = VT 
ONTIN 



1) VT = AINT(VT+.5) 



ih 



9 CONTINUE 

10 CONTINUE 

11 IF (NPR.LE.O) GO TO 12 
WRITE (6,51) I 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,I,FUN,CEN,I) 

12 lER =-l"’ 

GO TO 48 



90 



ono n on nno ooo o oonn o no no on 



c 

c 

c 



13 DO 14 J=1,NV 

SUMp) = §UM(J)+V(J,I) 

14 CEN(J) = SUMUj/FI 

TRY TO ASSURE FEASIBLE CENTROID FOR STARTING. 

NCE = NCE+1 

IF (KE(CEN) .EQ.O) GO TO 60 
SUmIj) = SUM(J) -V(J,I) 

GO TO 7 

60 NFE = NFE+1 

FUN (I) = FE(V(1,I)) 

15 CONTINUE 

END OF LOOP SETTING OF INITIAL COMPLEX. 

IF (NPR.LE.O) GO TO 17 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 

FIND THE WORST VERTEX, THE 'J’TH. 

J = 1 

DO 16 1=2, K 

IF_(FUN(J) .GE.FUN(I) ) GO TO 16 

16 CONTINUE 

BASIC LOOP. ELIMINATE EACH WORST VERTEX IN TURN, 
it must become NO LONGER WORST, NOT MERELY IMPROVED, 
find next- to-vertex, THE ' JN ' TH ONE. 

17 JN = 1 

IF (J.EQ.l) JN = 2 

DO 18 1 = 1, K 
IF (I.EQ.J) GO TO 18 
IF (FUNCJN) .GE.FUN(I) ) GO TO 18 
JN = I 

18 CONTINUE 

LIMT = NUMBER OF MOVES DURING THIS TRIAL TOWARD THE 
centroid DUE TO FUNCTION VALUE. 

LIMT = 1 

COMPUTE CENTROID AND OVER REFLECT WORST VERTEX. 



DO 19 1=1, NV 
VT = V(I,J) 

SUM(I) = SUM(I)-VT 

CENflJ = SUM(Ii/FKM 

VT = BETA''CEN?l)-ALPHA'-VT 

IF (IP.EQ.l) VT = AINT(VT+.5) 

INSURE THE EXPLICIT CONSTRAINTS ARE OBSERVED. 

19 V(I,J) = AMAX1(AMIN1(VT,BU(I) ) ,BL(I) ) 

NT = NT+1 

CHECK FOR IMPLICIT CONSTRAINT VIOLATION. 

20 DO 25 N=1,NLIM 
NCE = NCE+1 

IF (KE(V(1, J) ) .EQ.O) GO TO 26 

EVERY 'KV’TH TIME, OVER- REFLECT THE OFFENDING VERTEX 
through the BEST VERTEX. 

IF tMOD(N,KV) .NE.O) GO TO 22 
CALL FBV (K,FUN,M) 

DO 21 1=1, NV 

VT = BETA'W(I,M)-ALPHA*V(I, J) 
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IF (IP.EQ.l) VT = AINT(VT+.5) 

21 v(i,j) = aMxi(amini(vt,bu(i) ) , 

GO TO 24 

CONSTRAINT VIOLATION: 



BL(D) 



MOVE NEW POINT TOWARD CENTROID. 

22 DO 23 I=1,NV 

VT = .5''(CEN(I)+V(I, J) ) 

IF (IP.EQ.l) VT = AINT(vT+.5) 

V(I,j5 = VT 

23 CONTINUE 

24 NT = NT+1 

25 CONTINUE 

lER = 1 

CANNOT GET FEASIBLE VERTEX BY MOVING TOWARD CENTROID, 

OR BY OVER-REFLECTING THRU THE BEST VERTEX. 

IF (NPR.LE.O) GO TO 42 
WRITE (6,52j NT,J 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN, J) 

GO TO 42 

FEASIBLE VERTEX FOUND, EVALUATE THE OBJECTIVE FUNCTION. 

26 NFE = NFE+1 
FUNTRY = FE(V(1,J)) 

TEST TO SEE IF FUNCTION VALUE HAS NOT CHANGED. 

AFO = ABS (FUNTRY- FUNOLD) 

AMX = AMAX 1(ABS(EP" FUNOLD ) ,EP) 

ACTIVATE THE FOLLOWING TWO STATEMENTS FOR DIAGNOSTIC 
purposes only. 

WRITE (6,99) 

AFO .AMX , FUNTRY .FUNOLD , FUN ( J ) , FUN ( JN ) , NTFS , N 
9^ FORMAT (li,I3,6Ei5.7'2lS) 

IF (AFO.GT.AMi) GO t 6 27 
NTFS = NTFS+1 
IF (NTFS .LT.NCT) GO TO 28 
lER = 0 

IF (NPR.LE.O) GO TO 42 
WRITE (6,53) K 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 

GO TO 42 

27 NTFS = 0 

IS THE NEW VERTEX NO LONGER WORST? 

28 IF ( FUNTRY.lt. FUN (JN)) GO TO 34 

TRIAL vertex is STILL WORST: ADJUST TOWARD CENTROID. 
EVERY 'KV’TH time, OVER- REFLECT THE OFFENDING VERTEX 



C 

C 

C 

C 



C 

C 



through the BEST ^ERTEX 
LIMT = LIMT+1 
IF (MOD(LIMT.KV) .NE.O) GO TO 30 
CALL FBV (K,FUN,M) 

DO 29 1=1, NV 

VT = BETA"WU,M)-ALPHA*V(I,J) 

IF (IP.EQ.l) ^T = AINT(VT+.5) 

29 V(I,J) = AMAX1(AMIN1(VT,BU(I) ) ,BL(I) ) 

GO TO 32 



30 DO 31 I=1,NV 

VT = .5''(CEN(I)+V(I,J) ) 

IF (IP.EQ.l) VT = AiNT(VT+. 
V(I,J) = VT 



5) 



92 



nn nn nonnnon nn nno nnnnoooo non on 



C 

C 

c 

C 



31 CONTINUE 

32 IF (LIMT.LT.NLIM) GO TO 33 

CANNOT MA.KE THE ' J ' TH VERTEX NO LONGER WORST BY 
displacing toward 

THE CENTROID OR BY OVER-REFLECTING THRU THE BEST VERTEX. 
lER = 2 

IF (NPR .LE. 0) GO TO 42 
WRITE (6.52) NT, J 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN, J) 

GO TO 42 

33 NT = NT+1 
GO TO 20 

SUCCESS: WE HAVE A REPLACEMENT FOR VERTEX J. 

34 FUnTj) = FUNTRY 
FUNOLD = FUNTRY 
NPT = NPT+1 

EVERY 100 'TH PERMISSIBLE TRIAL, RECOMPUTE CENTROID 
summation to AVOID CREEPING ERROR. 

IF (MOD(NPT,100) .NE.O) GO TO 37 

DO 36 1=1, NV 
SUM (I) = 0. 

DO 35 N=1,K 

35 SUM(I) = §UM(I)+V(I,N) 

CEN(I) = SUM(I)/FK 

36 CONTINUE 

LC = 0 
GO TO 39 

37 DO 38 1=1, NV 

38 SUM(I) = §UM(I)+V(I, J) 

LC = J 

39 IF (NPR.LE.O) GO TO 40 

IF (MOD (NPT, NPR) .NE.O) GO TO 40 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,LC) 

HAS THE MAX. NUMBER OF TRIALS BEEN REACHED WITHOUT 
convergence? iF NOT, GO TO NEW TRIAL. 

40 IF (NT.GE.NTA) GO TO 41 

NEXT-TO-WORST VERTEX NOW BECOMES WORST. 

J = JN 
GO TO 17 

41 lER =3 ^ X 

IF (NPR.GT.O) WRITE (6,54) 

COLLECTOR POINT FOR ALL ENDINGS. 

1) CANNOT DEVELOP FEASIBLE VERTEX. lER = 1 

2) CANNOT DEVELOP A NO-LONGER- WORST VERTEX. lER = 2 

3) FUNCTION VALUE UNCHANGED FOR K TRIALS. lER = 0 

4) LIMIT ON TRIALS REACHED. lER = 3 

5) CANNOT FIND FEASIBLE VERTEX AT START. lER = -1 

42 CONTINUE 

FIND BEST VERTEX. 

CALL FBV Tk.FUN.M) ,, 

IF (IER.GE.3) GO TO 44 

RESTART IF THIS SOLUTION IS SIGNIFICANTLY BETTER THAN 
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no 



c 



c 

c 

c 



c 



c 



c 



c 



c 

c 



the previous, OR IF THIS IS THE FIRST TRY. 

IF (NPR.LE.O) GO TO 43 ^ _ 

WRITE (6,55) (M,YMN,FUN(M) ) 

43 IF (FUNTmJ .GE.YMN j do TO 47 ^ 

IF UbS(FUN(M)-YMN) .LE.AMAX1(EP,EP''YMN) ) GO TO 47 

GIVE IT ANOTHER TRY UNLESS LIMIT ON TRIALS REACHED. 

44 YMN = FUN?mT 
FUN(l) = M(M) 

DO 45 1=1, NV 



CEN(I 
SUM (I 
45 V(I,1 



^(I,M 



DO 46 1=1, NVT 

46 XS(I) = V(I,M) 

IF (IER.LT.3) GO TO 6 

47 IF (NPR.LE.O) GO TO 48 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,V(1,M) , -1) 
WRITE (6,56) ^(f^) 

48 RETURN 

49 FORMAT (50H0INDEX AND DIRECTION OF OUTLYING 
Ivariable at starti5) 

50 FORMAT (50H0IMPLICIT CONSTRAINT VIOLATED AT 
Istart. dead end.) 

51 FORMAT ('OCANNOT FIND FEASIBLE ', 14 ,’ TH VERTEX OR 
Icentroid at start.') 

52 FORMAT (lOHOAT TRIAL I4,54H CANNOT FIND FEASIBLE 
Ivertex which is no 

2L0NGER W0RST,I4,15X, 'RESTART FROM BEST VERTEX.') 

53 FORMAT ( 40H0FUNCTION HAS BEEN ALMOST UNCHANGED 
Ifor i5 , 7h trails) 

54 FORMAT (27H0LIMIT ON TRIALS EXCEEDED. ) 

55 FORMAT C^OBEST VERTEX IS NO. ',13, 'OLD MIN WAS',E15.7, 
1 ' NEW MIN IS ' ,E15 .7} 

56 FORMAT ( ' OMIN OBJECTIVE FUNCTION IS ',E15.7) 

END 

SUBROUTINE FBV (K,FUN,M) 

DIMENSION FUN(50) 

M = 1 

DO 1 1=2, K 

IF (FUN(M) .LE.FUN(I) ) GO TO 1 
M = I 

1 CONTINUE 

RETURN 

END 

SUBROUTINE BOUT (NT , NPT , NFE , NCE ,NV , NVT , V , K , FN , C , IK ) 
DIMENSION V(50, 50), FN(60),C(2i) 

WRITE (6,4) NT,NPT,NFE>CE 

DO 1 1=1, K 

WRI'TE (6,5) FN(I) , (V(J,I) , J=1,NV) 

IF (NVT.f,E.NV) GO TO 1 
NVP = NV+1 

WRITE (6,6) (V(J,I) , J=NVP,NVT) 

1 CONTINUE 

IF (IK.NE.O) GO TO 2 

WRITE (6,7) (C(I) ,I=1,NV) 

RETURN 

2 IF (IK.GE.O) GO TO 3 
WRITE (6,8) (C(I) ,I=liNV) 

RETURN 



94 



3 WRITE (6,9) IK, (C(I) ,I=1,NV) 

RETURN 

4 FORMAT ('ONO. TOTAL TRIALS = ',I5,4X, 

1 no. feasible trails = ',i5,4x, 

2 ; NO. FUNCTION EVALUATIONS = ',I5,4X, 

3 no. constraint evaluations = ' ,i5/ 

4'0 FUNCTION VALUE '6X INDEPENDENT VARIABLES/ 

SdependENT OR IMPLICIT CONSTRAINTS') 

5 FORMAT (IH , E18 . 7 , 2X , 7E14 . 7 / ( 21X , 7E14 . 7 ) ) 

6 FORMAT Uix!7E14.)) 

7 FORMAT aOH6CENTROID llX , 7E14 . 7 / ( 21X . 7E14 . 7 ) ) 

8 FORMAT ('0 BEST VERTEX ' , 7X , 7E14 . 7 7 ( 1 IX , 7E14 . 7 ) ) 

Q -nnnvxA-r > ' QCENTROID LESS ^X' ,I2,2X,7E14.7/ (21X, 

] 



9 FORMAT 
17el4.7) 

END 

FUNCTION FE(X 
DIMENSION Xn 
COMMON TDIFF 
CALL PLANT (X) 

FE=TDIFF 
RETURN 
END 

FUNCTION KE(X 
DIMENSION X(3 
KE = 0 
RETURN 
END 

//GO. SYS IN DD “ 

i THIS CARD SHOULD BE USED ONLY REGULAR SEA CASE 
//GO.FT12F001 DD DISP=SHR,DSN=MSS . S2160 . A341 



] 
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